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Abstract 

This  thesis  investigates  several  linear  and  nonlinear  feedback  control  methods  for 
satellite  formation  reconfigurations  and  compares  them  to  a  near  optimal  open  loop, 
discrete-time,  impulsive  maneuver.  The  reconfigurations  are  done  in  terms  of  a  set  of 
relative  parameters  that  define  an  orbit  about  the  leader  satellite  (or  center  reference  po¬ 
sition  if  a  leader  satellite  does  not  exist  at  the  center  of  the  formation).  The  purpose  of 
the  study  is  two-fold,  to  compare  the  control  usage  of  continuous  feedback  control  meth¬ 
ods  versus  a  discrete  burn  method  and  to  determine  if  nonlinear  control  techniques  offer 
significant  improvement  over  more  conventional  linear  control  laws.  Linear  Quadratic  Reg¬ 
ulators  (LQR),  LQR  with  linearizing  feedback,  State  Dependent  Riccati  Equation  (SDRE) 
and  sliding  mode  controllers  are  considered.  Simulations  showed  that  reconfigurations  for 
small  relative  orbits  were  adequately  controlled  using  linear  techniques. 
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A  STUDY  OF  LINEAR  VS.  NONLINEAR  CONTROL  TECHNIQUES  FOR  THE 
RECONFIGURATION  OF  SATELLITE  FORMATIONS 


I.  Introduction 

The  use  of  satellite  formations  to  do  the  work  of  larger  more  costly  single  satellites  has 
become  a  topic  of  great  interest  both  for  the  Air  Force  and  NASA.  The  potential  benefits 
of  using  formations  are  numerous.  Smaller,  cheaper  satellites  lead  to  reduced  launch  and 
life  cycle  costs.  Individual  satellites  can  be  replaced  giving  the  formation  the  ability  to 
upgrade  without  a  manned  mission  to  the  satellite.  Similarly,  spare  satellites  can  replace 
formation  components  as  they  fail.  This  creates  added  reliability  and  lifetime  to  a  mission. 
From  a  military  standpoint,  a  formation  has  a  greater  survivability  than  a  single  satellite 
as  multiple  satellites  are  often  harder  to  destroy  than  a  single  one.  Most  pertinent  to  this 
thesis  is  the  ability  of  a  formation  to  reconfigure  its  relative  geometry  in  order  to  meet 
changes  in  operational  needs,  adapt  to  the  failure  of  a  component,  or  allow  for  multi¬ 
mission  capability.  For  example,  a  formation  of  satellites  performing  a  mapping  mission 
using  synthetic  aperture  radar  (SAR)  could  change  its  ground  resolution  by  changing 
its  relative  positions  (14).  Although  satellite  formations  are  a  relatively  new  concept, 
several  applications  have  been  offered.  The  Air  Force  is  exploring  formations  for  use  in 
surveillance,  passive  radiometry,  terrain  mapping,  navigation,  communications,  and  ground 
target  identification  (22) .  In  the  civilian  sector,  satellite  formations  have  been  proposed  as 
a  way  to  provide  users  with  remote  sensing  data  in  the  visible  and  near  infrared  spectral 
bands  (2)  and  to  achieve  high  precision  localization  of  beacon  positions  (8).  Throughout 
this  thesis  the  reference  point  by  which  relative  motion  is  measured  is  called  the  ”  leader 
satellite”  or  simply  ’’leader”.  Satellites  occupying  relative  orbits  about  the  leader  satellite 
are  referred  to  as  ’’follower  satellites”  or  ’’follower”. 

From  a  dynamics  viewpoint,  two  of  the  larger  obstacles  to  functioning  satellite  for¬ 
mations  are  the  ability  to  keep  each  satellite  in  precise  relative  orbits  (stationkeeping)  and, 
if  the  need  arises,  to  change  those  relative  orbits  (reconfiguration).  Both  of  these  problems 
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require  accurate  relative  position  and  velocity  data  and  control  methods  that  use  an  ac¬ 
ceptable  amount  of  fuel.  To  produce  relative  position  and  velocity  data,  a  formation  could 
utilize  a  sensor  on  board  the  leader  satellite  (i.e.  bouncing  photons  off  follower  satellites), 
or  use  a  differential  global  positioning  system  (GPS)  method.  The  filtering  of  this  raw 
data  for  use  in  a  control  architecture  is  still  the  focus  of  intense  research  and  will  not  be 
examined  in  this  thesis.  The  control  method  will  take  the  relative  position  and  velocity 
data  provided  by  the  estimation  technique,  compare  it  to  the  desired  relative  orbit,  and 
use  that  error  to  calculate  the  needed  thrust  inputs  for  stationkeeping  and  reconfigura¬ 
tion.  This  leads  to  the  question:  How  will  desired  orbits  of  follower  satellites  be  specified? 
Mission  planners  will  most  likely  want  to  think  about  cluster  configurations  in  terms  of 
relative  orbits  about  a  leader  as  opposed  to  individually  specified  inertial  orbits  about  the 
Earth.  In  this  work,  relative  orbits  will  be  defined  in  terms  of  a  set  of  unique  parameters 
that  are  discussed  later.  Whatever  the  desired  relative  orbit,  the  follower  satellite  must 
also  be  in  a  Keplarian  orbit  in  the  absence  of  control  forces  and  perturbations.  To  keep  a 
satellite  in  a  non-Keplarian  orbit  requires  continuous  use  of  fuel,  an  unacceptable  situation 
for  a  satellite.  It  is  therefore  important  to  find  relative  orbits  that  not  only  fulfill  desired 
mission  needs  but  also  follow  a  Keplarian  orbit.  Only  orbits  that  follow  both  a  closed 
path  inertial  orbit  about  the  Earth  and  a  closed  path  relative  orbit  about  the  leader  in 
the  absence  of  control  forces  and  perturbations  will  be  examined  here.  The  stationkeeping 
controller  maintaining  the  relative  orbit  will  most  likely  be  separate  from  the  reconfigura¬ 
tion  controller  and  has  been  investigated  by  other  researchers  [(4), (27)].  Therefore,  this 
thesis  will  concentrate  on  linear  and  nonlinear  reconfiguration  controllers  with  total  fuel 
usage  (proportional  to  AV)  and  settling  time  as  the  key  parameters  for  comparison. 

1.1  Background 

The  basic  concepts  behind  relative  satellite  motion  were  developed  in  the  1960’s 
by  W.H.  Clohessy  and  R.S.  Wiltshire  who  were  working  on  methods  of  rendezvous  for 
orbital  construction  (6).  Although  developed  specifically  for  rendezvous,  these  relative 
equations  of  motion  work  well  in  the  design  of  useful  relative  orbits.  The  nonlinear  form 
of  the  Clohessy- Whiltshire  equations  (henceforth  known  as  CW  equations)  assumes  that 
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the  leader  satellite  exists  in  a  circular  orbit.  A  linear  form  of  the  CW  equations  can  be 
derived  if  the  additional  assumption  is  made  that  the  distance  between  the  leader  and 
follower  satellite  is  small  in  comparison  to  their  orbital  radii.  These  are  not  limiting 
factors  to  most  applications  of  satellite  formations.  One  notable  exception  is  the  use  of 
highly  eccentric  orbits  commonly  used  to  spend  more  time  over  a  certain  hemisphere  for 
communications  purposes  (i.e.  Molniya  orbits  (14)).  This  type  of  orbit  obviously  violates 
the  first  assumption  and  invalidates  the  CW  equations. 

Although  the  framework  for  developing  satellite  formations  has  been  around  since 
the  1960’s,  only  recently  has  this  problem  gained  widespread  attention  within  the  research 
community.  Early  research  tackled  the  estimation  problem.  The  Air  Force  sponsored  work 
that  determined  the  ability  of  an  on-board  U-D  covariance  Kalman  filter  algorithm  based 
on  the  CW  equations  to  estimate  relative  position  to  within  25  meters  (about  one  quarter  of 
a  radar  wavelength)  for  a  ten  satellite  orbital  formation.  Results  showed  the  position  error 
well  within  the  accuracy  requirement  (28).  This  research  was  expanded  to  include  the  use 
of  additional  measurement  data  from  other  satellites  in  the  formation  during  the  update 
cycle  which  produced  a  higher  degree  of  accuracy  in  the  relative  position  estimation  (7). 
A  later  Air  Force  Institute  of  Technology  (AFIT)  thesis  treated  the  Earth  as  an  imperfect 
sphere  as  opposed  to  a  point  mass  and  thus  developed  a  better  truth  model  which  included 
the  J2  term  of  the  Earth’s  geopotential  (17).  The  estimator  used  the  same  multiple  satellite 
data  and  Kalman  filter  as  used  in  (7).  Independent  investigation  of  an  on-board  relative 
position  estimator,  using  a  least-squares  error  analysis  of  short  range  satellite-to-satellite 
tracking  between  formation  members  was  examined  in  (25). 

Air  Force  interest  in  satellite  formations  is  apparent  in  the  TechSat  21  program 
which  was  created  to  explore  the  basic  technologies  required  to  create  a  viable  satellite 
formation  mission  (15).  The  TechSat  21  program  supports  research  into  collaborative 
behavior  of  satellites,  micro-propulsion,  micro-satellite  design,  micro-electro-mechanical 
systems,  smart  mechanisms,  multifunctional  structures,  and  lightweight  solar  arrays  (22). 
Some  researchers  have  investigated  the  parameterization  of  the  linearized  version  of  the 
CW  equations  (more  commonly  known  as  Hill’s  equations)  and  used  this  parameterization 
to  characterize  the  shapes  of  the  relative  orbits  (31).  Hill’s  equations  have  also  been  used  to 


3 


design  linear  controllers  to  solve  the  stationkeeping  problem  (24).  The  Air  Force  Research 
Laboratory  (AFRL)  has  supported  research  that  more  elegantly  defines  a  perceptive  frame 
which  supports  decentralized  feedback  (each  satellite  controls  its  own  relative  position 
error)  based  on  real-time  sensor  data  (11).  This  concept  is  expanded  and  used  with  an  LQR 
controller  on  several  simple  reconfigurations  (23).  The  research  into  satellite  formations  has 
matured  to  the  point  that  realistic  designs  of  the  intersatellite  links  have  been  proposed. 
These  are  the  communication  links  by  which  follower  satellites  would  talk  to  the  leader 
satellite  and,  if  need  be,  each  other  (9).  More  recent  attempts  to  define  the  problem  have 
concentrated  on  finding  equations  of  motion  that  are  not  constrained  to  circular  leader 
orbits  (30).  Similar  attempts  have  been  made  towards  the  estimation  of  relative  position 
in  elliptical  orbits  using  Kalman  filters  (10).  Although  it  is  desirable  to  be  able  to  generalize 
the  satellite  formation  problem  to  handle  elliptical  (and  thus  more  nonlinear)  leader  orbits, 
there  are  a  myriad  of  applications  for  the  constrained  case  which  is,  of  course,  much  simpler 
to  solve.  This  thesis  will  therefore  concentrate  on  the  constrained  case.  One  final  note 
about  satellite  formations;  the  European  Space  Agency  (ESA)  attempted  to  launch  the 
’’Cluster”  series  of  four  identical  satellites  aboard  the  maiden  flight  of  the  Ariane  5  as 
part  of  the  Inter-Agency  Solar  Terrestrial  Physics  program  (IASTP).  These  satellites  were 
designed  to  take  three-dimensional  field  and  plasma  measurements  of  the  magnetosphere; 
however,  the  launch  ended  in  failure  with  the  four  satellites  destroyed  in  June  1996.  Plans 
have  been  underway  since  November  1996  to  build  three  new  satellites  to  add  to  the 
preexisting  ground  spare  (3).  The  new  ’’Cluster”  satellites  were  launched  by  two  Russian 
Soyuz  rockets  (two  satellites  per  rocket)  on  16  July  and  9  August  2000  from  the  Baikonur 
Cosmodrome  in  Kazakhstan.  Their  science  mission  commenced  on  7  Feburary  2001  after 
an  extensive  checkout  (1). 

1.2  Objectives 

This  thesis  will  expand  on  previous  works  in  two  ways.  First,  it  will  compare  recon¬ 
figuration  simulations  using  both  linear  and  nonlinear  control  techniques  as  applied  to  the 
constrained  case  in  which  the  leader  is  in  a  circular  orbit.  It  is  important  to  determine  if 
nonlinear  control  techniques  provide  enough  of  a  benefit  (i.e.  lower  fuel  usage  and  recon- 
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figuration  times  versus  linear  techniques)  to  make  their  use  attractive  on  an  operational 
satellite.  If  not,  a  case  for  the  simpler  linear  techniques  will  be  made.  These  continuous 
control  techniques  will  also  be  compared  against  a  near  optimal,  open  loop,  discrete-time, 
impulsive  maneuver  (OLDTIM)  which  will  be  described  in  detail  later.  Second,  this  the¬ 
sis  will  employ  the  use  of  relative  parameters  (from  the  parameterization  of  the  linear 
equations  of  motion)  to  define  the  relative  orbit.  This  is  important  from  an  operational 
standpoint,  as  users  of  satellite  formations  will  most  likely  want  to  think  in  terms  of  relative 

parameters  are 
leader  satellite. 


orbits  rather  than  Keplarian  (Earth  centered)  orbits.  Therefore,  relative 
used  to  define  the  initial  and  final  relative  orbits  of  the  follower  about  the 
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II.  Relative  Satellite  Dynamics 

The  relative  equations  of  motion  are  based  on  the  inertial  orbit  equation  (Derivation  is  in 
Appendix  A) 

d  =  — =rr-*+  vc  +  vP  (1) 

\d\ 3 

where  //  is  the  gravitational  parameter  (for  Earth  y  =  398601  ^3-).  d  is  the  position 
vector  of  the  satellite,  vc  is  a  vector  of  control  forces  per  unit  mass,  and  vp  is  a  vector 
of  perturbation  accelerations  (oblate  Earth  effects,  air  drag,  solar  pressure,  etc.).  For 
now,  assume  that  the  satellite  operates  in  a  perturbations  free  environment  ( vp  =  0).  In 
’’Terminal  Guidance  System  for  Satellite  Rendezvous”,  (6)  Clohessy  and  Wiltwhire  used 
the  inertial  orbit  equation  to  derive  the  equations  of  motion  of  a  follower  satellite  about  a 
leader  satellite.  These  equations  are  the  foundation  upon  which  the  control  laws  developed 
here  are  based;  they  are  set  up  as  follows:  assume  a  moving  reference  frame  centered  on  the 
leader  satellite  which  is  in  a  circular  orbit  about  the  Earth.  The  Radial,  Cross  Track,  Out 
of  Plane  ( RCO )  coordinate  frame  is  formed  by  the  position  vector  of  the  leader  satellite 
in  inertial  space  [R],  the  velocity  vector  [C]  (which  for  a  circular  orbit  will  always  be 
perpendicular  to  the  position  vector),  and  the  cross  product  of  R  and  C  which  points 
out  of  the  orbit  plane  [i 6 ]  (same  direction  as  the  angular  momentum  vector  of  the  leader 
satellite)  as  shown  in  Figure  1.  In  the  RCO  frame,  the  follower’s  relative  equations  of 


Figure  1  The  RCO  Frame 
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motion  are  [(6), (23)]  (Derivation  is  in  Appendix  C): 


r  —  2  uc  —  u2{k0  +  r) 


1  - 


k±. 


c  +  2  u>r  —  u2c 
o  +  u )20 


1  - 


[(k0+r)2+c2+o2]i  J 

-  vCc  =  o 


Vc  =  0 


[(fco+r)2+c2+o2]^  j 


(2) 


L[(fco+r)2+c2+o2]5  J 


-  vCo  =  0 


where  r  is  the  relative  position  in  the  radial  direction,  c  is  the  relative  position  in  the  cross 
track  direction,  and  o  is  out  of  the  leader  orbit  plane.  The  remaining  two  variables  are 
properties  of  the  leader  orbit:  u>  is  the  angular  frequency  (or  mean  motion)  of  the  leader 
orbit  and  k0  is  the  radius  of  the  leader  orbit  where  the  relationship  between  the  two  is 
based  on  Kepler’s  third  law  (29): 

ui  = 

Equations  2  are  used  in  the  development  of  nonlinear  controllers,  but  can  be  linearized  to 
facilitate  the  use  of  linear  control  methods.  When  the  distance  between  the  follower  and 
leader  is  small  compared  to  the  radii  of  their  orbits,  linearization  yeilds  (Derivation  is  in 
Appendix  D): 


r  —  2  uic  —  3  u2r  —  uCr  =  0 

c  +  2ur  -  vCc  =  0  (4) 

o  +  uj2o  -  vCo  =  0 

This  is  the  more  familiar  form  of  the  Clohessy-Wiltshire  results  commonly  known  as  Hill’s 
equations.  If  the  control  vector  vc  is  set  to  zero,  the  homogeneous  form  of  Hill’s  equations 
can  be  parameterized  as  (31)  (Derivation  is  in  Appendix  E): 

r  =  p  sin  (ut  +  9)  +  a 

c  =  2pcos(wt  +  6)  —  +  b  (5) 

o  =  mpsm(u)t  +  6)  +  2np  cos(wf  +  6) 
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where 


a  —  2c0+4a;rn 

U) 

(6) 

u  _  ojc0— 2rp 

UJ 

(7) 

=  \J(r0-a)2+(%f 

(8) 

">=*$ 

(9) 

_  _  o0r0u-\-6oUj{a-ro) 

J  ~~  2[r^+o;:i(a— r0)2] 

(10) 

6  =  arctan  [w(r;-a)] 

(11) 

which  fall  out  of  the  derivation  of  Equations  5  (Note  that  the  subscript  o  denotes  initial 
conditions).  In  the  simulations  run  for  this  thesis,  a,  fr,  and  p  are  defined  in  kilometers,  9 
in  degrees,  and  m  and  n  in  terms  of  a  unitless  slope  (Appendix  E.l).  What  we  first  notice 
about  the  parameterization  is  that  the  second  equation  has  a  term  that  grows  linearly  with 
time.  Since  u  will  never  be  zero,  it  is  easy  to  see  that  a  must  always  equal  zero  or  else  the 
relative  orbit  will  not  follow  a  closed  path.  Thus  there  is  an  initial  condition  constraint 
that  c0  =  —2 u>r0.  Further  analysis  of  the  follower  orbit  can  be  found  in  Section  2.1.  Taking 
the  derivatives  of  Equations  5  to  find  the  relative  velocity  yields: 

r  —  puj  cos(u )t  +  6) 

c  —  —2  pu  sin  {ujt  +  9)  —  --a  (12) 

b  —  mpu  cos(c Jt  +  9)  —  2  npu  sin  (cvt  +  9) 

Equations  5  and  12  can  be  used  to  convert  desired  relative  orbits  (in  terms  of  relative 
parameters)  into  inertial  position  and  velocity  at  a  given  time.  First  the  RCO  frame 

relationship  to  the  inertial  frame  ( UK )  is  needed  to  find  the  transformation  matrix 
(j  Fr  ameA—2— Frame  B 


IJK Vector  =  cRC°2IJK  *RCO  Vector 
CRC02IJK  =  [R\C\0] 


(13) 

(14) 
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where  R,  C,  and  6  are  the  unit  vectors  for  each  direction  expressed  in  the  UK  frame.  Note 
that  in  the  UK  system,  /  is  in  the  direction  of  the  first  point  of  Aires,  J  is  perpendicular 
to  I  in  the  Earth’s  equatorial  plane  and  K  =  I  x  J  (29).  If  the  leader  satellite’s  inertial 
position  and  velocity  are  IJKL  and  IJKL  respectively  and  the  follower’s  inertial  position 
and  velocity  are  IJK M  and  IJKM,  all  expressed  in  the  UK  frame,  then  the  radial  unit 
direction  is  simply  the  unit  vector  of  the  leader’s  position 


0  =  RxC 


(17) 


To  find  the  inertial  position  of  the  follower  in  the  UK  frame,  the  relative  position  vector 
is  transformed  to  UK  and  added  to  the  leader’s  position 

r 

ijkM  =  cRC°2IJK  c 

o 

To  find  the  inertial  velocity  we  need  to  take  the  inertial  derivative  of  the  position  vector. 
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thus  the  inertial  velocity  of  the  follower  satellite  becomes 


IJK$l  _  QRC02IJK 


r  —  cw 
c  +  rcu 


o 


(25) 


Going  from  the  inertial  to  the  relative  frame  is  a  simple  manipulation  of  the  above  equations 
and  is  shown  in  Appendix  F. 


2.1  The  Relative  Orbit 

A  manipulation  of  Equations  5  yields  better  insight  into  what  path  the  follower  will 
take  about  the  leader  and  what  each  parameter  represents.  If  a  is  set  to  zero,  squaring  the 
first  two  equations,  adding  them  together,  and  using  the  trigonometric  identity  sin(r)2  + 
cos(t)2  =  1  yields 

r2  +  —  — -  =  p2[sin(u;t  +  9)]2  +  p2[cos(u>t  +  0)]2  =  p2  (26) 


Dividing  both  sides  by  p2 


Substituting  r  =  psm(cvt  +  9)  and  (c  — 


(c-6)^ 
+  4p2 


b )  =  2pcos(u)t  +  0)  into  the  third  equation 


(27) 


o  =  mr  +  n(c  —  b ) 


(28) 


Equation  27  shows  that  the  path  of  the  follower  satellite  in  the  radial/cross-track  plane  will 
be  an  ellipse  with  a  semi  minor  axis  equal  to  p  and  a  semi  major  axis  equal  to  2 p.  Thus  p 
defines  the  size  the  relative  orbit.  The  b  term  defines  an  offset  in  the  cross-track  direction 
and  shows  that  a  closed  path  relative  orbit  can  exist  without  the  leader  satellite  located  at 
the  center  of  the  ellipse.  In  fact,  the  center  of  the  relative  orbit  can  exist  anywhere  along 
the  velocity  direction  of  the  leader  provided  it  is  not  so  far  from  the  leader  as  to  render 
Hill’s  equations  invalid.  The  m  and  n  terms  in  Equation  28  define  the  slope  of  the  orbit 
of  the  radial/out-of-plane  and  cross- track/out-of-plane  planes  respectively.  The  9  term  is 
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the  initial  angle  with  respect  to  the  cross  track  direction  and  t  is  time.  Figure  2  gives  a 
graphical  representation  of  Equation  27.  The  two  foci  of  the  ellipse  fall  on  the  major  axis 


Figure  2  Path  of  Follower  in  R/C  Plane 


a  distance  /  from  the  center  of  the  ellipse  where  /  is  found  via  (13) 


/2  =  (Major  Axis)2  —  (Minor  Axis)2 


(29) 


Thus 

f  =  =  =  (30) 

The  eccentricity  of  an  ellipse  is  found  with  the  following  relationship  (13). 

e=-rJ-  ,  =^  =  4  (31) 

Major  Axis  2 p  2 

Thus  the  eccentricity  of  the  elliptical  path  in  the  R/C  plane  will  always  be  If  m  is  set 
to  zero,  the  follower  track  in  the  C/O  plane  will  be  a  line  of  slope  n.  Likewise,  if  n  is  set 
to  zero,  the  follower  track  in  the  R/O  plane  will  be  a  line  of  slope  m. 

Other  than  the  constraints  on  a  (to  prevent  relative  orbit  drift),  6,  and  p  (to  maintain 
the  validity  of  the  linear  equations)  a  satellite  formation  mission  planner  is  free  to  chose  any 
set  of  parameters  needed  to  create  the  relative  orbit  required  by  the  mission.  By  choosing 
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a  unique  set  of  relative  parameters,  the  user  is  choosing  a  unique  set  of  initial  conditions. 
The  tradeoff  for  the  ease  of  relative  parameters  is  that  the  command  signal  is  now  based 
on  the  linear  relative  equations  of  motion,  thus  even  at  steady  state  conditions  (i.e.  no 
reconfiguration),  there  will  always  be  error  between  the  inertially  propagated  system  and 
the  commanded  input.  This  error  is  quantified  and  discussed  in  the  following  section. 

2.2  Error  Using  Linear  vs.  Nonlinear  Equations 

Linearizing  the  CW  equations  requires  the  assumption  that  the  magnitude  of  the 
relative  position  vector  is  small  compared  to  the  radius  of  the  leader  orbit.  In  order  to 
make  valid  conclusions  about  reconfigurations  using  a  linear  based  command  signal,  it  is 
important  to  quantify  the  error  produced  by  using  the  linear  equations  as  opposed  to  the 
nonlinear.  The  simulation  shown  in  Figure  3  is  used  to  quantify  that  error.  The  simulation 


Figure  3  Error  Analysis  Simulation  Setup 

was  run  for  various  values  of  p  (which  defines  the  size  of  the  relative  orbit)  and  leader  orbit 
radii.  Relative  parameters  were  used  to  determine  initial  conditions  for  both  the  linear  and 
nonlinear  equations  of  motion  (EOMs)  via  Equations  5  and  12.  These  initial  conditions 
are  then  used  in  both  sets  of  equations  (2  and  4)  and  propagated  through  two  complete 
relative  orbits.  The  error  is  the  absolute  value  of  the  difference  in  the  magnitude  of  the 
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position  vectors.  The  relative  parameters  for  each  simulation  are  shown  in  Table  1.  There 


Sim  1 

Sim  2 

Sim  3 

Sim  4 

Sim  5 

Sim  8 

Sim  9 

Sim  10 

p  (km) 

0.1 

0.2 

0.3 

0.5 

0.8 

1 

1.4 

1.8 

2.4 

3 

a  (km) 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

9  (deg) 

45 

45 

45 

45 

45 

45 

45 

45 

45 

45 

b  (km) 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

m  (unitless) 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

n  (unitless) 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

Table  1  Absolute  Error  Simulation  Relative  Parameters 


were  three  leader  orbit  radii  used;  a  low  Earth  orbit  (semi  major  axis  =  6800  km),  a 
geosynchronus  orbit  (semi  major  axis  =  42241  km),  and  an  orbit  in  between  (semi  major 
axis  =  17720  km).  For  each  orbit  radius,  p  was  varied  from  0.1  to  3  km.  The  R  and  O 
directions  showed  periodic  error  patterns  as  shown  in  Figures  4  and  5  respectively  where 
each  line  corresponds  to  a  specific  p  and  the  error  increases  as  p  increases.  The  error  in 

Radius  of  Leader  Orbit  =  6800  (km) 


Figure  4  R  Direction  Absolute  Position  Error 


the  C  direction  has  a  secular  growth  element  and  is  shown  in  Figure  6.  Note  that  the  cross 
track  error  is  an  order  of  magnitude  greater  than  in  the  other  two  directions.  The  presence 
of  this  secular  growth  indicates  that  setting  a  =  0  (via  c0  =  —  2u>r0),  does  not  eliminate  all 
secular  terms  in  the  cross  track  direction  (at  least  for  the  nonlinear  equations).  It  is  this 
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Time  (min) 


Figure  5  O  Direction  Absolute  Position  Error 


Radius  of  Leader  Orbit  =  6800  (km) 


Figure  6  C  Direction  Absolute  Position  Error 
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growth  in  the  error  of  the  C  direction  that  is  the  greatest  concern.  For  one  thing,  creating 
initial  conditions  based  solely  on  the  relative  parameters  will  not  guarantee  a  stable  relative 
orbit  (that  is  an  orbit  that  does  not  drift  over  time).  From  a  classical  orbital  elements 
standpoint,  this  means  that  the  follower  and  leader  satellites  have  a  different  semi-major 
axis  and  thus  different  periods  translating  to  a  drift  in  their  relative  position.  This  is 
much  more  of  a  factor  for  stationkeeping  controllers  than  reconfigurations  since  the  former 
occurs  over  long  time  periods  thus  allowing  cross  track  error  to  build  up.  Depending 
on  the  settling  time  of  the  controller,  the  cross  track  error  during  reconfiguration  may 
be  acceptable.  Looking  at  the  maximum  error  (peak  of  each  absolute  error  curve)  as  a 
function  of  p  produces  Figures  7,  8,  and  9.  As  expected,  the  larger  the  leader  satellite  orbit 


Figure  7  R  Direction  Max  Error 

and/or  p,  the  smaller  the  max  error  between  the  linear  and  nonlinear  equations.  Note  that 
in  the  R  and  O  directions,  the  max  error  graph  is  valid  for  any  number  of  relative  orbits 
but  in  the  C  direction  the  graph  would  be  shifted  upwards  for  subsequent  revolutions  due 
to  the  growth  term.  In  any  case,  the  C  direction  error  will  be  the  determining  factor  in 
deciding  how  large  p  can  be. 
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Figure  8  C  Direction  Max  Error 


Figure  9  0  Direction  Max  Error 


2.3  Simulation  Outline 


The  reconfiguration  simulation  used  to  obtain  results  for  this  thesis  is  setup  according 
to  Figure  10.  Equations  5  and  12  are  used  to  convert  relative  parameters  and  current  time 


Figure  10  Reconfiguration  Simulation  Diagram 


into  a  relative  position  and  velocity  (in  the  RCO  frame).  This  is  the  command  by  which 
the  relative  measurement  data  is  subtracted  in  order  to  produce  an  error  signal.  It  is 
assumed  that  relative  position  and  velocity  data  is  available  (i.e.  full  state  feedback).  The 
error  signal  is  multiplied  by  the  controller  gain  to  find  the  control  inputs.  These  inputs 
(vc)  are  converted  to  the  I JK  frame  and  used  in  the  inertial  orbit  equation 


+  Vc  +  Vp 


(32) 


where  d  is  the  position  of  the  satellite  with  respect  to  the  Earth  in  the  inertial  frame  and  vp 
is  the  perturbation  acceleration.  The  orbits  of  both  the  leader  and  follower  are  propagated 
via  Equation  32  (note  that  for  the  leader  vc  =  0)  to  find  a  new  inertial  position  and  velocity 
vector  for  the  leader  and  follower  which  are  differenced  and  transformed  to  the  RCO  frame 
(Appendix  F  explains  the  process).  This  is  the  measurement  data  compared  against  the 
command  closing  the  loop.  Although  it  may  not  be  realistic  to  expect  real  time  position  and 
velocity  data,  the  purpose  of  this  research  is  to  compare  controller  designs,  not  estimators. 
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Future  work  will  include  estimators  that  effectively  deal  with  noise  and  propagate  states 
between  measurements.  The  results  from  each  controller  will  also  be  compared  against 
the  near  optimal  open  loop,  discrete-time  impulsive  maneuver  (OLDTIM)  (Section  3.1). 
No  claim  is  made  that  the  discrete  burn  method  used  is  optimal,  only  that  it  represents  a 
recognized,  efficient  maneuver  for  orbit  transfer. 

The  inertial  orbit  equation  is  propagated  via  the  Dormand-Prince  pair  (an  explicit 
Runge-Kutta  formula).  Since  both  follower  and  leader  orbits  are  either  circular  or  near 
circular,  there  is  little  need  for  a  variable  step  solver.  A  fixed  time  step  method  was  chosen 
to  give  the  simulation  user  more  control  over  the  run  time  of  the  simulation.  A  time  step 
of  15  seconds  was  selected  for  this  thesis  based  on  analysis  of  the  integration  error  as  a 
function  of  time  step  and  because  it  provided  reasonable  run  times.  The  simulation  used 
to  propagate  the  leader  and  follower  inertial  states  was  initialized  for  a  circular  orbit  and 
run  for  one  revolution.  The  position  error  is  the  maximum  error  between  the  magnitude 
of  the  position  vector  and  the  semi  major  axis  of  the  circular  orbit  (used  to  find  initial 
conditions)  and  is  shown  in  Figure  11.  The  error  simulation  shows  that  for  a  15  second 


Figure  11  Integration  Error  vs.  Time  Step  (Position) 

time  step,  the  error  in  the  magnitude  of  the  position  vector  is  under  1  x  10-6.  The  velocity 
error  is  the  maximum  error  between  the  magnitude  of  the  velocity  vector  and  the  expected 
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speed  of  the  satellite  based  on  the  equation 


V=M  (33) 

V  fco 

where  k0  is  the  orbital  radius  (semi  major  axis).  The  velocity  error  is  shown  in  Figure  12. 
The  error  simulation  shows  that  for  a  15  second  time  step,  the  error  in  the  magnitude  of 


Time  Step  (s) 


Figure  12  Integration  Error  vs.  Time  Step  (Velocity) 
the  velocity  vector  is  under  2  x  10-9. 

The  rest  of  the  Simulink  simulation  (Appendices  K,  L,  M,  and  N)  involves  converting 
the  data  into  relative  parameters,  classical  orbital  elements,  and  error  signals. 
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III.  Reconfiguration  Control  Laws 

This  chapter  describes  each  of  the  four  control  techniques  used  to  develop  controllers  for 
the  relative  orbit  reconfiguration  problem.  The  control  techniques  used  for  this  thesis 
are  the  linear  quadratic  regulator  (LQR),  a  nonlinear  LQR  technique  using  linearizing 
feedback,  state  dependent  riccati  equations  (SDRE),  and  sliding  mode  control.  The  first 
section  explains  the  assumptions  and  derives  the  equations  used  to  calculate  the  OLDTIM 
maneuver  AV  and  settling  time.  This  is  the  open  loop  maneuver  that  will  be  used  as  a 
comparison  point  for  the  continuous  feedback  methods. 

3.1  The  Open  Loop,  Discrete-Time,  Impulsive  Maneuver 

The  open  loop,  discrete-time,  impulsive  maneuver  (OLDTIM)  consists  of  four  sepa¬ 
rate  burns.  Experience  with  several  simulation  runs  shows  that  all  of  the  follower  satellites’ 
classical  orbital  elements  (COEs)  change  during  the  reconfiguration  except  for  the  semi 
major  axis.  As  discussed  before,  changes  in  the  semi  major  axis  would  cause  the  relative 
orbit  to  drift  as  it  would  affect  the  period  of  the  follower  satellites’  inertial  orbit.  True 
anomaly  is  not  considered  as  it  describes  the  position  of  a  satellite  within  an  orbit,  not  the 
size,  shape,  or  attitude  of  the  orbit.  The  change  in  eccentricity  can  be  performed  using  a 
Hohmann  transfer  which  is  an  optimal  two  burn  maneuver  for  a  circular  to  circular  orbit 
transfer  and  near  optimal  for  eccentric  orbits.  The  inclination  and  longitude  of  the  ascend¬ 
ing  node  can  be  changed  with  a  single  burn  given  certain  assumptions  presented  later.  The 
final  parameter,  argument  of  perigee  is  changed  with  a  single  burn  as  well.  The  author 
is  aware  that  other  more  optimal  methods  for  discrete  burns  are  available,  however,  this 
method  allows  for  the  fewest  assumptions  and  is  significantly  easier  to  calculate.  Given 
the  fact  that  the  OLDTIM  maneuver  is  only  a  comparison  point  and  not  meant  to  bound 
an  optimal  transfer,  the  tradeoff  is  acceptable. 

3.1.1  The  Hohmann  Maneuver  (Eccentricity  Change).  The  Hohmann  maneuver 
will  perform  the  eccentricity  change.  To  calculate  the  necessary  AV,  consider  Figure  13 
where  ra  is  the  radius  of  apogee  of  the  larger  orbit  and  rp  is  the  radius  of  perigee  of  the 
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Figure  13  Hohmann  Maneuver  Setup 

smaller  orbit.  These  values  are  related  to  the  semi  major  axis  (A)  and  eccentricity  (e)  (29) 


Ai  =  A2  =  A  (34) 

I'd  =  A2(l  +  e2)  =  ^4(1  +  e2)  (35) 

rp  =  Ai{l  -  ei)  =  A(1  -  ex)  (36) 

The  total  AV  for  the  Hohmann  maneuver  will  be 

AV Hohmann  =  AVx  +  A  Vy  (37) 


To  calculate  the  difference  in  orbital  velocity  at  points  X  and  Y,  we  need  to  find  the  velocity 
of  the  satellite  at  each  point  in  its  orbit  (29) 


V  = 


/ 


(38) 
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3.1.2  The  Inclination/LAN  Maneuver  (Plane  Change).  The  second  maneuver 
will  change  the  inclination  and  longitude  of  the  ascending  node  through  a  plane  change. 
Consider  Figure  14  in  which  we  assume  that  both  the  initial  and  final  orbits  are  circular 
with  the  same  radius.  This  is  a  good  assumption  as  long  as  the  relative  orbit  is  small  thus 
requiring  a  small  eccentricity  for  the  follower  satellite.  All  possible  orbits  for  a  certain 
orbital  radius  exist  on  the  surface  of  a  sphere.  In  Figure  14,  \I>  represents  the  angle 
between  the  planes  of  the  two  orbits,  f 2  is  the  longitude  of  the  ascending  node,  and  i  is  the 
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Figure  14  Inclination/Longitude  of  the  Ascending  Node  Maneuver  Setup 
inclination.  With  spherical  trigonometry  (5) 

=  arccos[cos(ii)  cos(z2)  +  sin(ii)  sin(i2)  cos(fi2  —  fli)]  (44) 

With  the  angle  between  the  two  planes  known,  the  AV  can  be  calculated.  Consider  Figure 
15  where  Vq  is  the  speed  of  the  satellite  (constant  for  all  points  on  the  surface  of  the 
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is  found  via  (5) 


Vo  =  \lv 


\2_ 

k„ 


ko 


(45) 


where  ka  is  the  radius  of  the  circular  orbit.  Using  the  trigonometric  formula  for  the  sine 
of  an  angle  (Figure  15) 

✓  «tv  \  1  at; 

(46) 


sin(!)= 


\^Vlnc[LAN_ 

Vo 


Solving  for  AV 

£VVInc/LAN  =  2Vo  sin  ^  =  2y^-sin  j  (47) 

Since  this  is  a  single  instantaneous  burn  and  will  happen  between  the  initial  and  final 
burns  of  the  Hohmann,  the  settling  time  is  zero. 


3.1.3  The  AP  Maneuver  (Orientation  Change).  The  final  discrete  burn  will 
change  the  angle  of  the  line  of  perigee  and  is  set  up  in  Figure  16.  Since  only  the  argument 


Figure  16  AP  Maneuver  Setup 
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of  perigee  (AP)  changes,  all  other  orbital  elements  will  remain  the  same,  thus 


ei  =  e2  =  e  (48) 

Ai  =  A2  =  A  (49) 

Since  the  semilatus  rectum  (p)  is  a  function  of  e  and  A 

Pi=P2  =  P  =  A(1-  e2)  (50) 

Zooming  in  on  Figure  16  the  following  relationships  can  be  seen  (Figure  17) 


Figure  17  AP  Maneuver  Setup  Zoom 


a  =  Pi 


(51) 


V 1  =  27T  -* 


A  AP 
2 


(52) 


^2 


A  AP 
2 


(53) 


where  (3  is  the  angle  between  the  position  and  velocity  vectors  of  each  orbit  and  v  is  the 
angle  between  the  line  of  perigee  and  the  position  vector  (also  known  as  true  anomaly). 
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Taking  the  cosine  of  both  sides  of  Equations  52  and  53  (26) 


,  .  /  A AP\  ,n  ,  / AAP\ 

cos(i/i)  =  cos  I  2tt - - — I  =  cos(27r)  cos  I — - — j  +  sin(27r) 

/A  AP\ 

=  COS  [—) 

cos(i/2)  =  cos  j 


Thus 


cos  ( i/i )  =  cos  (1/2) 
AAP 

z/i  =  i/2  =  1/  =  — — 


Finally,  looking  at  how  (3  is  constructed  in  Figure  18  we  can  say  that 


Figure  18  f3  Relationships 

A  =  f  -  71 
A  =  f  +  72 


where  7  is  the  flight  path  angle  and  is  defined  as 


7  =  arccos 


where  Vjv  is  the  velocity  normal  to  the  position  vector  (Vr  is  the  velocity  parallel  to  the 
position  vector).  Substituting  Equations  58  and  59  into  Equation  51 

( 7 r 

«  =  (  2 +  T2 


)  -  (|  -  7i)  =  72  +  7i 


(61) 


All  the  angle  relationships  are  now  in  place  to  solve  for  the  AV  of  the  AP  change.  Note 
on  Figure  16  that  at  the  maneuver  point,  the  two  orbits  share  a  common  position  vector. 
Relating  the  position  vector  to  conic  section  geometry  (29) 


\di\  =  \d2\  =  d  = 


A(  1  -  e2) 

1  +  e  cos(^) 


P 

1  +  e  cos(i/) 


(62) 


Equation  38  shows  that  the  magnitude  of  the  velocity  is  a  function  of  the  magnitude  of 
the  position  vector  and  the  semi  major  axis,  both  of  which  are  identical  for  orbit  one  and 
two.  Solving  for  V  and  substituting  in  Equations  50  and  62 


\Vi\  =  \v2\  =  v  = 


1  +  2e  cos(i')  +  e2 
P 


(63) 


Talcing  the  cross  product  of  the  two  position  vectors  with  their  respective  velocity  vectors 
yields  the  angular  momentum 


\di  x  Vr | 


\d2xV2\  = 


\Hi\  =  dV sin(/?i)  =  dV sin  -  71^ 

=  dV  sin  ^0  cos (71)  -  cos  ^0  sin(7i) 
\H2\  =  dV  sin (/?2)  =  dV  sin  +  72^ 

in  ^0  cos(72)  +  cos  ^0  sin(72) 


=  dV 


The  magnitude  angular  momentum  ( H )  is  (5) 

\H\  =  VMP 


dV  cos  (71) 

(64) 

dV  cos(72) 

(65) 

(66) 
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Since  p\  =  P2 


II 

teli 

to 

(67) 

dV  cos  (71)  =  dV  cos  (72) 

(68) 

i- 

II 

f5 

II 

i-H 

(69) 

By  Equation  60  this  means  that 

V/Vl  =  Vn2  =  Vn 

(70) 

Substituting  7  into  Equation  61 

a  =  27 

(71) 

and  relating  7  to  angular  momentum 

dV  cos(7)  =  \H\  =  y/jlp 

(72) 

cos(7)  =  W 

(73) 

Substituting  Equation  60  for  cos(7) 

II 

a* 

(74) 

Cancelling  the  V  and  substituting  Equation  62  for  d 

%  =  VW1  +  eCOSM=^[  1  +  ecosM] 

(75) 

Taking  the  cosine  of  Equation  71  and  employing  the  double  angle  cosine  trigonometric 

identity  (26) 

cos(a)  =  cos(27)  =  2[cos(7)]2  —  1 

(76) 

Substituting  Equation  60  for  7 


cos(o;)  =  2 


( 

\VN 1 

'll 

2 

r^vi 

cos  1  arccos 

v  J 

). 

-1  =  2 

v  J 

-  1  = 


2VN2  -  V2 
V2 


(77) 
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Now  that  we  have  c*,  Vi,  and  V2  the  law  of  cosines  states  that 


A  VAp2  =  V?  +  Vi  -  2V1V2  cos(a)  =  2V2(1  -  cos  (a))  (78) 

Substituting  for  cos(a)  (Equation  76) 

AVap2  =  2V2  1  ~  2VNy^  ~  =W2-Wn2  (79) 

Expanding  this  with  the  values  of  V  and  Vjv  via  Equations  63  and  75 

1  +  2e  cos(i/)  +  e2  a.  ,l2  /0(V. 

- - qi  +  ecos  i/  2  (80) 

P  P  J 

Adding  the  fractions,  expanding  the  squared  term,  and  simplifying 

AW2  =  4„  [ At  ~  =  4',e2lsin(l/)l2  (81) 

P  J  P 

Finally,  talcing  the  square  root  of  both  sides  and  substituting  AML  for  v  (per  Equation  57) 

AVAP  =  2e^sin(^)  (82) 

or  in  terms  of  classical  orbital  elements 

AVap  =  2e{: i(i^)sin  (83) 

3.1.4  Total  AV  and  Settling  Time.  Adding  Equations  42,  47,  and  83  to  get  the 
total  AV  for  the  OLDTIM 

A  VoLDTIM  =  AVnohmann  +  AVInc/LAN  +  AV^p  (84) 

To  determine  the  total  maneuver  time  (or  settling  time)  visualize  the  following  sequence. 
The  first  burn  will  place  the  spacecraft  on  the  Hohmann  transfer  orbit  and  occurs  at 
apogee  of  the  initial  orbit.  There  are  two  opportunities  to  do  the  Inclination/LAN  burn 
located  180  degrees  apart,  therefore  one  opportunity  will  occur  between  apogee  and  perigee 


Ay4P2  =  4  /i 
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and  will  not  affect  the  maneuver  time.  The  third  burn  occurs  at  perigee  and  places  the 
spacecraft  onto  the  final  orbit.  The  fourth  and  final  burn  changes  the  argument  of  perigee. 
Figure  16  shows  that  the  spacecraft  either  travels  an  angle  of  from  the  third  burn  at 
perigee  till  the  fourth  burn  or  the  fourth  burn  occurs  before  the  final  Hohmann  burn.  In 
the  form  case,  the  following  relationship  calculates  the  travel  time  (29) 


STap  —  At  — 


i—[E-esm(E)) 

jji 


where  E  is  the  eccentric  anomaly  and  is  equal  to  (29) 


(85) 


E  =  2  arctan 


where  v  =  Thus  the  ’’settling  time”  of  the  OLDTIM  maneuver  is 


(86) 


STqlDTIM  =  ST  Hohmann  +  STap 


(87) 


STqldtim  = 


(88) 


3.2  Linear  Quadratic  Regulator 

The  Linear  Quadratic  Regulator  (LQR)  is  based  on  linear,  time-invarient  systems  of 
the  form 

x  =  Ax  +  Bu  (89) 

where  re  is  a  vector  of  states  and  u  is  a  vector  of  control  inputs.  LQR  minimizes  the 
quadratic  performance  index 

J  =  i  J  ( xtQx  +  uTRuj  dt  (90) 
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where  Q  is  the  state  weighting  matrix  and  R  is  the  control  weighting  matrix.  For  the 
relative  orbit  reconfiguration  problem,  the  states  represent  the  relative  position  and  velocity 
in  the  RCO  frame  and  thus  all  states  are  weighed  equally  giving  Q  a  diagonal  form  in  which 
each  diagonal  element  is  the  same  magnitude.  Likewise,  excess  control  usage  is  equally 
bad  in  all  three  directions  thus  R  has  a  diagonal  form  in  which  each  diagonal  element  is 
the  same  magnitude.  High  values  of  the  elements  in  the  state  weighting  matrix  result  in 
faster  movement  from  initial  to  desired  states.  High  values  of  the  elements  in  the  control 
weighting  matrix  result  in  lower  control  usage.  For  the  purposes  of  this  simulation,  Q  is 
set  to  the  identity  matrix  and  the  elements  of  R  will  be  varied.  Placing  Hill’s  equations 
(4)  in  the  form  of  Equation  89  yields 
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With  A,  B,  Q,  and  R  an  algebraic  Riccati  equation  is  solved  to  find  S 


'Vf 

Vc 


Vo 


(91) 


SA  +  ATS  +  Q-  SBR~1BtS  =  0 


(92) 


The  static  gain  controller  matrix  is  then 


K  =  R~1BtS 


(93) 
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The  control  inputs  will  therefore  be  (18) 


where  the  subscript  com  represents  the  commanded  states. 


(94) 


3.3  Linear  Quadratic  Regulator  With  Linearizing  Feedback 

In  order  to  use  a  linear  technique  such  as  LQR  without  losing  the  nonlinear  aspects 
of  the  dynamics,  we  can  use  linearizing  feedback  (LF).  This  technique  is  useful  only  if  the 
nonlinear  equations  can  be  split  into  linear  and  nonlinear  parts.  Standard  linear  techniques 
can  then  be  applied  to  that  part  of  the  problem  while  nonlinear  terms  are  added  in  after 
control  is  calculated.  Starting  with  the  nonlinear  CW  equations  (2),  assume  the  control 
inputs  ( vTc ,  vCci  v0c )  have  the  form 


r  —  2uc  —  3c o2r  —  v'T  =  0 

c  +  2  ujf  —  v'c  =  0  (96) 

d  +  u>2o  —  vfo  —  0 


33 


The  static  gain  controller  matrix  is  found  using  the  same  LQR  process  as  outlined  in 
Section  3.2  but  now  the  linearizing  feedback  must  be  added  to  the  output  of  the  controller 
to  find  the  control  signal 


" 

VCr 

Vcc 
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.  V°°  . 

-u2(k0  +  r ) 


1  - 


JL 


—uj2c 


[(fco+r,)2+c2+o2]^  J 


+  3c o2r 


1- 


— c o2o 


[(k0+r)2+c2  +o2\2 


*L 


L  [{ko+r)2+c2+02]^  J 
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r  r  com 
T  T com 
C  —  Ccom 
C  ~  ccom 
°  Ocom 
0  —  Ocom 


(97) 


where  the  subscript  com  represents  the  commanded  states. 


3-4  State  Dependent  Riccati  Equations 

The  State  Dependent  Riccati  Equation  (SDRE)  technique  uses  the  nonlinear  CW 
equations. 


r  —  2  ojc  +  uj2(k0  +  r) 


1  - 


£ 


[(ko+r)2+c2+02]^  J 


c  =  —2  cvf  +  u2c 
0  =  — uj2  o 

Assume  the  following  simplification 

(Jq  — 
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Gf 
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crc 


Substituting  these  equations  into  Equations  98  yields 

r  =  2ujc  +  u2arr  +  v ^ 
c  =  -2c or  +  u)2<jcc  +  vCc 


(98) 


(99) 

(100) 

(101) 


(102) 
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6  =  -tJ2a0o  +  vCo 


In  matrix  form 
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(103) 


Since  oy,  <rc,  and  cr0  change  at  each  time  step,  the  SDRE  technique  continuously  updates 
Equation  103  and  finds  a  new  controller  gain  matrix  via  the  Riccati  equation  discussed  in 
Section  3.2  and  uses  the  new  gain  to  find  the  control  signal. 


r  -  room 


VCr 

Vcc 


=  -K 


Vc0 


r  ~~  T com 
c  Ccom 
C  ~~  ccom 
0  °com 


0  °com 


(104) 


where  the  subscript  com  represents  the  commanded  states. 

For  LQR,  LQR  with  linearizing  feedback,  and  SDRE,  the  control  weight  ( R )  will  be 
varied  for  each  simulation  run.  Trial  runs  show  that  a  range  of  about  109  to  1013  produce 
acceptable  results.  Lower  than  109  results  in  unacceptably  high  control  usage  while  larger 
than  1013  translates  to  unacceptable  settling  times.  The  specifics  can  be  found  in  Chapter 
IV. 


3.5  Sliding  Mode  Control 

Sliding  mode  control  (SMC)  is  a  nonlinear  control  technique  for  imprecise  systems. 
Since  the  commanded  input  is  based  on  linear  equations  (creating  an  imprecise  system) 
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this  type  of  controller  is  a  candidate.  SMC  is  derived  by  defining  a  time  varying  surface 
upon  which  the  system  has  some  desired  dynamics  (21)  and  is  defined  as 

/  d  \n_1 

s(x,  t)  =  +  Aj  x  (105) 

where  x  is  the  tracking  error 

F  T com 
T  ~  T com 
C  “  ccom 
c  Ccom 
0  ~  °com 
°  —  °com 

where  the  subscript  com  represents  the  commanded  states,  n  is  the  order  of  the  system 
(in  this  case  n  —  2)  and  A  are  the  real  poles  that  determine  the  system  dynamics  on  the 
surface  s(x,t).  Since  the  system  order  in  this  case  is  two,  Equation  105  reduces  to 

s  =  x  +  Xx  (107) 

The  goal  of  SMC  is  to  drive  s  to  zero 

s  —  x  +  Xx  =  0  (108) 

The  sliding  mode  function  from  the  Nonlinear  Synthesis  toolbox  from  Optimal  Synthesis 
Inc  was  used  to  implement  sliding  mode  control  for  this  simulation.  The  function  required 
a  pole  (A)  for  each  direction,  a  parameter  77  which  designates  the  rates  at  which  the  system 
will  converge  to  the  sliding  surface,  and  e  which  designates  the  boundary  layer  of  the  sliding 
surface  (16).  Within  the  boundary  layer  defined  by  e,  no  control  is  used. 

Trial  runs  using  sliding  mode  control  on  the  nonlinear  CW  equations  showed  the  best 
rj  and  e  are 

rj  =  10~5'8 
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c  =  10~6 


Instead  of  varying  the  control  weight  ( R )  as  in  LQR  and  SDRE  techniques,  the  pole 
locations  are  varied  for  each  sliding  mode  run.  Trial  runs  showed  that  the  poles  should 
exist  no  farther  from  the  ju>  axis  than  -0.003  for  acceptable  control  usage  and  settling  time 
results. 

The  Nonlinear  Synthesis  toolbox  requires  a  formulation  of  the  problem  such  that  the 
states  are  driven  to  zero;  therefore  the  states  need  to  be  modified  so  that  they  are  equal 
to  the  error  signal  instead  of  the  relative  position  and  velocity  of  the  follower  satellite. 
Starting  with  the  nonlinear  CW  equations 
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The  modified  states  and  their  derivative  are 
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Making  the  same  simplifications -as  in  the  SDRE  derivation 
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yields 


r  =  2  ujc  +  u>2arr  +  v ^ 

c  =  —  2ur  +  u2acc  +  vCc  (1 14) 

o  =  -u2a0o  +  vCo 

Substituting  the  commanded  and  actual  states  in  Equations  109  (note  that  vCcom 

T’ com  ?  =  [2 0JCcom  4"  w  crrComrcom]  [2cUC  +  U)  <J rT  +  VCr] 
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—  —  2u)(rcom  —  r)  +  uj  (o’Ccomccom  crcc)  —  vCc  (116) 
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(115) 


Ccom  0 


[~V2cr ocomOCom\  ~  [-W2cr00  +  V0c] 
-V2(crocomOcom  ~  VqO)  ~  VQc 


(117) 


The  Nonlinear  Synthesis  toolbox  requires  a  function  file  that  passes  the  derivative  of  the 
state  ( x )  given  the  current  state  (#)  and  a  vector  of  controls  (vc),  thus  this  function  would 
pass  back  the  following  vector 
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IV.  Results 


The  following  chapter  presents  a  comparison  of  the  different  control  techniques  for  sev¬ 
eral  different  types  of  reconfigurations.  Important  information  about  how  well  each  con¬ 
troller  performs  the  reconfiguration  can  also  be  found  by  examining  the  control  usage  and 
state/COE  time  histories  as  well  as  other  data  produced  during  each  simulation.  For 
this  reason,  each  control  technique  has  an  appendix  where  the  complex  parameter  change 
reconfiguration  (Section  4.4)  is  presented  and  discussed  in  more  detail. 

•  LQR  in  Appendix  G 

•  LQR  with  linearizing  feedback  in  Appendix  H 

•  SDRE  in  Appendix  I 

•  Sliding  Mode  in  Appendix  J 

The  settling  time  for  each  maneuver  is  measured  from  the  start  of  the  reconfiguration 
not  t  =  0,  likewise,  the  AV  listed  for  each  reconfiguration  is  the  total  AV  measured  at 
the  settling  time  minus  the  total  AV  measured  at  the  reconfiguration  start  time.  This 
means  the  AV  does  not  include  any  of  the  steady  state  control  energy  used  to  correct  the 
error  between  the  linear  and  nonlinear  equations  before  or  after  the  reconfiguration.  Of 
course  due  to  the  linear  command  signal  there  will  always  be  additional  error  that  the 
controller  is  forced  to  deal  with  as  shown  in  Section  2.2.  The  first  section  of  this  chapter 
will  characterize  this  error  and  discuss  its  effects  on  AV  for  each  controller. 

4.1  The  Impact  of  Using  Linear  vs.  Nonlinear  Equations  on  AV 

In  the  reconfiguration  simulation,  the  command  signal  is  based  on  the  parameterized 
Hill’s  equations  (linear  EOMs)  while  the  orbits  are  propagated  inertially  (nonlinear  EOMs). 
This  causes  an  error  between  the  commanded  and  the  actual  relative  positions  even  without 
a  reconfiguration  and  requires  control  usage.  It  is  this  control  usage  (and  subsequent  AV 
costs)  that  is  of  particular  interest  and  must  be  quantified.  Running  the  reconfiguration 
simulation  without  a  change  in  the  commanded  relative  parameters  is  used  to  examine 
the  AV  costs  of  the  indicated  errors.  Using  the  relative  parameters  listed  in  Table  2,  a 
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simulation  was  run  using  a  time  step  of  15  seconds.  The  slope  for  each  compensator  is 


P  (km) 

0.5 

a  (km) 

0 

6  (deg) 

45 

b  (km) 

0 

ra  (unitless) 

1 

n  (unitless) 

0 

Table  2  AV  Error  Simulation  Relative  Parameters 

determined  via  a  least  squares  line  fit  and  is  listed  in  Table  3.  The  use  of  AV  to  null  out 
the  command  error  is  fairly  linear  for  all  values  of  control  weight.  For  the  LQR  and  SDRE 
techniques,  smaller  values  of  R  (which  translate  to  higher  control  usage)  resulted  in  higher 
slopes  while  the  reverse  was  true  for  LQR  with  LF  (Figure  19).  The  converse  was  true  for 
larger  values  of  R  (Figure  20).  Experience  shows  that  the  total  AV  expended  during  a 


I 

> 

< 


Figure  19  Gain  Methods  Steady  State  AV  Costs  R  =  109 

reconfiguration  using  LQR,  LQR  with  LF,  or  SDRE  is  several  orders  of  magnitude  greater 
than  the  AV  used  to  correct  for  the  linear  command.  For  this  reason  and  the  fact  that 
these  three  compensators  have  similar  slopes,  the  excess  AF  usage  will  be  accepted.  The 
sliding  mode  steady  state  costs  are  much  more  significant.  Using  the  parameters  in  Table 
2,  the  slopes  for  sliding  mode  are  listed  in  Table  4.  It  is  obvious  that  the  sliding  mode 


Time  (days) 
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Table  3  Steady  State  AV  Slopes  (Gain  Methods) 


controller  is  using  a  large  amount  of  control  energy  during  steady  state  conditions.  This 
is  due  to  excess  chatter  on  the  sliding  mode  surface;  a  known  problem  with  sliding  mode 
control.  The  negative  effects  of  this  chatter  will  become  obvious  in  the  example  simulations 
and  is  the  main  reason  that  the  sliding  mode  technique  requires  higher  A  Vs  than  LQR, 
LQR  with  LF,  and  SDRE.  The  slopes  seem  to  increase  as  the  sliding  mode  pole  moves 
away  from  the  ju  axis.  Figures  21  and  22  show  the  steady  state  AV  cost  for  the  highest 
(P=-0.003)  and  lowest  (P=-0.0006)  SMC  slopes  respectively. 
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re  21  Pole  Method  Steady  State  AV^  Costs  (P  =  —0.003 


Time  History  of  Delta  V  for  Poie=-0.0006  [Rho=0.5  (km)] 


Figure  22  Pole  Method  Steady  State  AT^  Costs  (P 
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AV  Usage  (™  per  day) 

LQR 

Pole=-0.0006 

358.9385 

Pole=-0.0012 

361.1704 

Pole=-0.0018 

363.4057 

Pole=-0.0024 

364.9747 

Pole=-0.003 

365.577 

Table  4  Steady  State  AV  Slopes  (Pole  Method) 

4-2  Simple  p  Reconfiguration  (Small  Relative  Orbit) 

This  reconfiguration  changed  only  the  parameter  p  thereby  changing  the  semi  ma- 
jor/minor  axes  of  the  relative  orbit.  The  choice  of  the  initial  and  final  p  are  small  enough 
that  the  linear  Hill’s  equations  are  very  close  to  the  nonlinear  CW  equations.  From  the 
error  analysis  done  in  Section  2.2  ,  we  expect  no  more  than  a  1  meter  error  in  the  R  and 
O  directions  and  up  to  10  meters  in  the  C  direction  (depending  on  settling  time)  for  a 
10,000  km  leader  orbit.  The  initial  and  final  relative  parameters  are  given  in  Table  5. 
The  simulation  uses  a  15  second  time  step  with  the  leader  in  a  circular  orbit  with  classi¬ 
cal  orbital  elements  (COEs)  listed  in  Table  6.  The  AV  results  for  controllers  in  which 


Initial 

Final 

p  (km) 

0.5 

1.5 

a  (km) 

0 

0 

e  (deg) 

45 

45 

b  (km) 

0 

0 

m  (slope) 

1 

1 

n  (slope) 

0 

0 

Table  5  Small  Orbit  p  Simulation  Initial  and  Final  Parameters 


Semi-Major  Axis 

10000  km 

Eccentricity 

0 

Inclination 

10  deg 

Argument  of  Perigee 

0  deg 

Longitude  of  the  Ascending  Node 

0  deg 

Initial  True  Anomaly 

10  deg 

Table  6  Small  Orbit  p  Simulation  Leader  COEs 


control  weight  was  varied  are  listed  in  Table  7.  The  results  show  that  for  higher  values 
of  R,  the  two  nonlinear  methods  do  provide  a  lower  AV  as  compared  to  LQR,  however, 
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AV  (m/s) 

OLDTIM 

LQR 

LQR  With  LF 

SDRE 

13.8110215 

13.8110284 

13.8110049 

1 

7.8518765 

7.8519592 

7.8518429 

4.2552826 

4.2552369 

4.2551942 

2.3860346 

2.3855585 

2.3859874 

R  =  1013 

1.7847727 

1.7842614 

1.7847499 

Table  7  Small  Orbit  p  Simulation  AV  Results  (Control  Weight  Methods) 

this  improvement  is  on  the  order  of  10-6  and  thus  negligible.  The  AV  results  for  slid¬ 
ing  mode  control  (which  varied  pole  placement)  are  listed  in  Table  8.  The  sliding  mode 


AV  (m/s) 

OLDTIM 

SMC 

Pole=-0.0006 

0.632092 

36.5888133 

Pole=-0.0012 

0.632092 

19.9482038 

Pole=-0.0018 

0.632092 

15.3333261 

Pole=-0.0024 

0.632092 

14.2883201 

Pole=-0.003 

0.632092 

13.8828796 

Table  8  Small  Orbit  p  Simulation  AV  Results  (Pole  Method) 

method  shows  a  significant  increase  in  control  energy  as  opposed  to  the  gain  techniques. 
Section  4.1  shows  that  compared  to  the  gain  techniques,  sliding  mode  control  expends  a 
large  amount  of  control  energy  even  during  steady  state  conditions  due  to  the  chattering 
effect.  The  controller  is  required  to  expand  even  more  energy  on  top  of  this  to  do  the 
reconfiguration  and  thus  at  a  disadvantage  compared  to  the  gain  techniques. 

The  OLDTIM  AV  can  be  broken  out  into  its  three  separate  discrete  burn  maneuvers 
and  is  shown  in  Table  9.  The  changes  in  the  COEs  of  the  follower  satellite  upon  which 
the  OLDTIM  maneuver  is  based  are  listed  in  Table  10.  The  OLDTIM  break  out  shows 


Hohmann  AV 
Inc/LAN  AV 
AP  AV 

0.000417502  (m/s) 
0.631353  (m/s) 
0.000320607  (m/s) 

Total 

0.632092  (m/s) 

Table  9  Small  Orbit  p  Simulation  OLDTIM  AV 

that  a  majority  of  the  control  energy  in  the  discrete  burn  is  used  to  change  the  inclination 
and  longitude  of  the  ascending  node  even  though  they  are  relatively  small  angle  changes. 
Plane  changes  are  usually  more  AV  intensive  than  coplanar  maneuvers. 
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Initial 

Final 

A 

Semi-Major  Axis  (km) 

10000 

10000 

N/A 

Eccentricity 

0.0000499951 

0.000149973 

0.0000999777 

Inclination  (deg) 

10.0023 

10.007 

0.00469583 

Argument  of  Perigee  (deg) 

304.987 

304.929 

0.0581968 

Longitude  of  the  Ascending  Node  (deg) 

-0.00945911 

-0.0283561 

0.018897 

Table  10  Small  Orbit  p  Simulation  Follower  COEs  Changes 


The  settling  time  results  are  listed  in  Tables  11  and  12.  The  criteria  for  settling  time 
is  that  the  position  error  for  each  direction  in  the  RCO  frame  falls  and  remains  within 
+/-  10  meters.  It  is  assumed  that  at  this  point  a  stationkeeping  controller  will  take  over 
and  reduce  or  maintain  the  error  to  required  levels.  As  expected  after  seeing  the  AV 


Settling  Time  (min) 

OLDTIM 

LQR 

LQR  With  LF 

SDRE 

82.933410 

20.00 

20.00 

20.00 

82.933410 

36.50 

36.50 

36.50 

82.933410 

62.00 

62.00 

62.00 

82.933410 

107.50 

107.25 

107.50 

82.933410 

287.75 

288.00 

287.75 

Table  11  Small  Orbit  p  Simulation  Settling  Time  Results  (Control  Weight  Methods) 


Settling  Time  (min) 

OLDTIM 

SMC 

Pole=-0.0006 

82.933410 

146.00 

Pole=-0.0012 

82.933410 

82.25 

Pole=-0.0018 

82.933410 

65.25 

Pole=-0.0024 

82.933410 

65.00 

Pole=-0.003 

82.933410 

69.00 

Table  12  Small  Orbit  p  Simulation  Settling  Time  Results  (Pole  Method) 

costs,  the  linear  and  nonlinear  gain  controllers  continue  to  perform  virtually  identically. 
To  see  a  few  other  trends,  this  data  is  presented  graphically  in  Figures  23  and  24  for  AV 
and  Figures  25  and  26  for  settling  time.  Notice  that  for  the  control  weight  methods,  the 
total  AV  for  the  reconfiguration  approaches  the  OLDTIM  AV  nearly  asymptotically  as 
R  increases.  Also  note  that  as  AV  decreases,  settling  time  increases  (a  normal  tradeoff  in 
designing  controllers).  This  is  not  the  case  for  sliding  mode  control  where  AF  and  settling 
time  decrease  as  the  poles  move  from  the  origin  to  around  -0.002  on  the  real  axis.  Since  the 
steady  state  AV  usage  of  sliding  mode  is  so  high,  the  quicker  it  settles  to  the  commanded 
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orbit,  the  less  AV  it  requires.  The  large  amount  of  energy  expended  by  the  sliding  mode 
controller  during  steady  state  is  overwhelming  its  ability  to  control  the  follower  satellite 
at  a  reasonable  AV  cost.  To  compare  all  the  methods  together  and  determine  their 


Gain  Methods 


Figure  23  Small  Orbit  p  Simulation  AV  (Control  Weight  Methods) 


Pole  Methods 


Figure  24  Small  Orbit  p  Simulation  AV  (Pole  Method) 

best  operating  point  ( R  or  pole  location),  each  method  is  graphed  via  AV  versus  settling 
time  (Figure  27) .  Note  that  the  closer  the  curve  comes  to  the  origin,  the  better  the  overall 
performance  of  the  method.  Once  a  specific  method  is  determined,  a  designer  can  find  the 
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Gain  Methods 


Figure  25 


Small  Orbit  p  Simulation  Settling  Time  (Control  Weight  Methods) 


Pole  Methods 


Figure  26  Small  Orbit  p  Simulation  Settling  Time  (Pole  Method) 


Method  Comparison 


Figure  27  Small  Orbit  p  Simulation  Method  Comparison 

point  at  which  the  tradeoff  between  AV  and  settling  time  is  acceptable  based  on  mission 
parameters  (on  board  fuel,  acceptable  off  line  times,  frequency  of  reconfigurations,  etc.). 

For  this  reconfiguration  in  which  p  is  changed  and  the  relative  orbit  is  small  compared 
to  the  radius  of  the  leader  orbit,  LQR,  LQR  with  linearizing  feedback,  and  SDRE  work 
significantly  better  than  sliding  mode  control  (again,  mostly  due  to  chatter  problems). 
Each  of  these  three  methods  produces  virtually  identical  results,  thus  there  is  no  advantage 
to  using  nonlinear  techniques  over  linear  ones  in  this  case.  This  is  not  a  surprising  result 
as  the  relative  orbit  is  small  enough  that  the  linear  equations  of  motion  are  a  very  good 
approximation. 

4-3  Simple  p  Reconfiguration  (Large  Relative  Orbit) 

The  next  relative  orbit  reconfiguration  changed  the  parameter  p  by  the  same  amount 
as  in  the  first  simulation,  however,  in  this  case  the  initial  and  final  p  are  much  larger 
thus  creating  more  error  between  the  linear  and  nonlinear  EOMs.  For  this  larger  p,  error 
analysis  shows  that  we  can  expect  an  almost  1200  meter  error  in  the  R  and  O  directions 
and  up  to  14,000  meters  in  the  C  direction  (depending  on  settling  time)  for  a  10,000  km 
leader  orbit.  The  initial  and  final  relative  parameters  are  given  in  Table  13.  The  simulation 
uses  a  15  second  time  step  with  the  leader  in  a  circular  orbit  with  classical  orbital  elements 
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(COEs)  listed  in  Table  14.  The  AV  results  for  controllers  in  which  control  weight  was 


Initial 

Final 

p  (km) 

40.5 

41.5 

a  (km) 

0 

0 

6  (deg) 

45 

45 

b  (km) 

0 

0 

m  (slope) 

1 

1 

n  (slope) 

0 

0 

Table  13  Large  Orbit  p  Simulation  Initial  and  Final  Parameters 


Semi-Major  Axis 

10000  km 

Eccentricity 

0 

Inclination 

10  deg 

Argument  of  Perigee 

0  deg 

Longitude  of  the  Ascending  Node 

0  deg 

Initial  True  Anomaly 

10  deg 

Table  14  Large  Orbit  p  Simulation  Leader  COEs 


varied  are  listed  in  Table  15.  The  control  usage  in  this  simulation  is  much  different  than  in 


AP  (m/s) 

OLDTIM 

LQR 

LQR  With  LF 

SDRE 

R  =  10y 

1.260372 

15.6921806 

13.6182920 

15.6912682 

R  =  1010 

1.260372 

11.4366934 

7.7630169 

11.4385323 

R  =  1011 

1.260372 

9.3192450 

4.7136968 

9.3163821 

R  =  1012 

1.260372 

9.0168279 

3.9794505 

9.0305245 

R  =  1013 

1.260372 

10.1362775 

8.1898422 

10.1407846 

Table  15  Large  Orbit  p  Simulation  AV  Results  (Control  Weight  Methods) 

the  near  linear  case  of  the  first  simulation.  The  LQR  and  SDRE  controllers  showed  near 
identical  A  Vs,  but  LQR  with  LF  proved  to  use  significantly  less  total  control  energy.  The 
AF  results  for  sliding  mode  control  (which  varied  pole  placement)  are  listed  in  Table  16. 
As  before,  sliding  mode  control  shows  a  significant  increase  in  expended  control  energy  as 
compared  to  the  gain  techniques  due  to  the  chattering  effect  discussed  in  Section  4.1. 

The  AV  for  each  of  the  OLDTIM  maneuvers  is  shown  in  Table  17.  The  changes  in 
the  COEs  of  the  follower  satellite  upon  which  the  OLDTIM  maneuver  is  based  are  listed 
in  Table  18.  The  OLDTIM  breakout  is  significantly  different  than  the  first  simulation. 
Both  the  Hohmann  and  AP  maneuvers  require  two  orders  of  magnitude  more  control 
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AV  (m/s) 

OLDTIM 

SMC 

Pole=-0.0006 

1.260372 

38.8888040 

Pole=-0.0012 

1.260372 

19.6637421 

Pole=-0.0018 

1.260372 

16.1290137 

Pole=-0.0024 

1.260372 

14.7930510 

Pole=-0.003 

1.260372 

14.3854720 

Table  16  Large  Orbit  p  Simulation  AV"  Results  (Pole  Method) 


Hohmann  AV 
Inc/LAN  Al/ 
AP  AV 

0.182102  (m/s) 
0.667614  (m/s) 
0.410657  (m/s) 

Total 

1.260372  (m/s) 

Table  17  Large  Orbit  p  Simulation  OLDTIM  AV" 

energy  while  the  plane  change  AV  only  marginally  increased.  Note  that  their  is  more 
error  associated  with  the  initial  and  final  COEs  as  they  too  are  a  calculated  from  the 
linear  command  signal.  Thus  the  OLDTIM  AV"  calculation  has  a  larger  margin  of  error 
as  well. 

The  settling  time  results  are  listed  in  Tables  19  and  20.  The  criteria  for  settling  time 
is  that  the  position  error  for  each  direction  in  the  RCO  frame  falls  and  remains  within 
+/-  10  meters  after  which  a  stationkeeping  controller  will  take  over  and  reduce  the  error 
to  required  levels.  An  analysis  of  the  position  error  history  shows  that  only  LQR  with 
LF  and  sliding  mode  control  are  able  to  reduce  the  position  error  to  the  prescribed  +/-  10 
meters.  The  error  plots  of  LQR  and  SDRE  show  that  the  error  will  oscillate  with  a  fixed 
amplitude  (Figure  28)  about  the  zero  axis  with  no  apparent  dampening.  The  settling  times 
of  LQR  and  SDRE  are  consequently  approximately  equal  to  the  simulation  run  time  and 
meaningless.  Analysis  of  the  position  error  at  different  gains  indicates  that  the  amplitude  is 
a  function  of  the  control  weight  R .  Lower  values  of  R  (less  restrictive  on  control  usage)  will 
allow  both  controllers  to  keep  error  within  +/-  10  meters  but  at  unacceptable  AV"  costs. 
To  see  the  trends  in  AV  and  settling  time,  this  data  is  presented  graphically  in  Figures  29 
and  30  for  AV  and  Figures  31  and  32  for  settling  time.  The  inability  of  LQR  and  SDRE 
to  damp  the  position  and  velocity  state  error  means  that  AV"  will  not  approach  OLDTIM 
value.  LQR  with  LF  also  does  not  approach  OLDTIM  asymptotically  and  instead  appears 
to  have  a  local  minimum  at  around  R  =  1012.  As  before,  the  settling  time  increases  as  R 
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Initial 

Final 

A 

Semi-Major  Axis  (km) 

10000 

10000 

N/A 

Eccentricity 

0.00401814 

0.0041296 

0.000111458 

Inclination  (deg) 

10.1916 

10.1968 

0.00517635 

Argument  of  Perigee  (deg) 

303.981 

303.054 

0.927511 

Longitude  of  the  Ascending  Node  (deg) 

-0.743455 

-0.761246 

0.0177914 

Table  18  Large  Orbit  p  Simulation  OLDTIM  AV 


Settling  Time  (min) 

OLDTIM 

LQR 

LQR  With  LF 

SDRE 

■HSiHi 

82.946683 

90.50 

20.00 

90.50 

R  =  1010 

82.946683 

158.25 

36.50 

158.25 

82.946683 

216.00 

62.00 

216.00 

R  =  1012 

82.946683 

331.00 

107.25 

331.00 

82.946683 

576.00 

288.00 

576.00 

Table  19  Large  Orbit  p  Simulation  Settling  Time  Results  (Control  Weight  Methods) 


Settling  Time  (min) 

OLDTIM 

SMC 

Pole=-0.0006 

82.946683 

152.25 

Pole=-0.0012 

82.946683 

79.75 

Pole=-0.0018 

82.946683 

68.50 

Pole=-0.0024 

82.946683 

66.75 

Pole=-0.003 

82.946683 

70.75 

Table  20  Large  Orbit  p  Simulation  Settling  Time  Results  (Pole  Method) 


LQR  R=1e+013 


Figure  28  LQR  Position  Error  for  the  Large  Orbit  p  Simulation  (R  =  1013) 
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increases.  Sliding  mode  also  damps  out  the  error  signal  to  zero  but  at  a  AV  cost  an  order 
of  magnitude  above  LQR  with  LF  and  appears  to  asymptotically  approach  a  minimum 
AV  value  as  the  pole  moves  from  the  origin.  Unlike  LQR  with  LF,  the  AV  decreases 
with  shorter  settling  times  due  to  the  high  steady  state  AV  cost  discussed  in  Section  4.1. 

Graphing  all  of  the  methods  together  for  comparison  and  optimal  point  determination 


Gain  Methods 


Figure  29  Large  Orbit  p  Simulation  AV  (Control  Weight  Methods) 


Pole  Methods 


Figure  30  Large  Orbit  p  Simulation  AV  (Pole  Method) 
yields  Figure  33.  Note  that  the  closer  the  curve  comes  to  the  origin,  the  better  the  overall 
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Gain  Methods 


R  (Control  Weight) 


Figure  31  Large  Orbit  p  Simulation  Settling  Time  (Control  Weight  Methods) 


Pole  Methods 


Figure  32  Large  Orbit  p  Simulation  Settling  Time  (Pole  Method) 
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Method  Comparison 


Figure  33  Large  Orbit  p  Simulation  Method  Comparison 
performance  of  the  method. 

Although  the  change  in  relative  parameters  for  this  simulation  is  identical  to  the  first 
simulation,  it  is  clear  that  using  linear  techniques  will  not  be  adequate.  The  command 
signal  is  based  on  the  linear  EOMs  and  thus  the  expended  control  energy  will  increase  as 
the  relative  orbit  size  becomes  a  significant  percentage  of  the  orbit  radius.  Only  LQR  with 
LF  and  sliding  mode  were  able  to  perform  the  reconfiguration  to  within  the  specified  error 
tolerance  and  of  these  two,  LQR  with  LF  did  the  maneuver  with  significantly  less  AV  and 
in  an  acceptable  settling  time. 

4-4  Complex  Parameter  Reconfiguration 

This  simulation  presents  a  more  general  relative  orbit  reconfiguration  in  which  all 
parameters  change  except  for  a  and  b  (thus  leaving  the  leader  satellite  at  the  center  of  the 
formation) .  Like  the  first  simulation,  the  choice  of  initial  and  final  p  is  small  enough  that 
the  linear  Hill’s  equations  are  a  very  close  approximation  of  the  nonlinear  CW  equations. 
The  small  relative  position  magnitude  means  that  we  expect  no  more  than  a  1  meter  error 
in  the  R  and  O  directions  and  10  meters  in  the  C  direction  for  a  10,000  km  orbit.  The 
initial  and  final  relative  parameters  are  given  in  Table  21.  The  simulation  uses  a  15  second 
time  step  with  the  leader  in  a  circular  orbit  with  classical  orbital  elements  (COEs)  listed 
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in  Table  22.  The  AV  results  for  controllers  in  which  control  weight  was  varied  are  listed 


Initial 

Final 

P  (km) 

0.5 

1.5 

a  (km) 

0 

0 

e  (deg) 

45 

60 

b  (km) 

0 

0 

m  (slope) 

1 

1.5 

n  (slope) 

0 

1 

Table  21  Complex  Simulation  Initial  and  Final  Parameters 


Semi-Major  Axis 

10000  km 

Eccentricity 

0 

Inclination 

10  deg 

Argument  of  Perigee 

0  deg 

Longitude  of  the  Ascending  Node 

0  deg 

Initial  True  Anomaly 

10  deg 

Table  22  Complex  Simulation  Leader  COEs 


in  Table  23.  The  results  show  once  again  that  the  gain  techniques  are  nearly  equal  in 


AV  (m/s) 

OLDTIM 

LQR 

LQR  With  LF 

SDRE 

R=  109 

2.351894 

31.4307824 

31.4307616 

31.4307437 

R  =  1010 

2.351894 

17.2212144 

17.2213265 

17.2211453 

R  =  1011 

2.351894 

9.2636521 

9.2637386 

9.2635397 

R  =  1012 

2.351894 

5.2854682 

5.2845003 

5.2854414 

R  =  1013 

2.351894 

4.2184381 

4.2181246 

4.2185121 

Table  23  Complex  Simulation  AV  Results  (Control  Weight  Methods) 


their  reconfiguration  performance  and  no  one  method  stands  out.  The  AV^  results  for 
sliding  mode  control  (which  varied  pole  placement)  are  listed  in  Table  24.  The  sliding 
mode  method  once  again  requires  a  much  greater  AV  expenditure  as  compared  to  the 
gain  methods  due  to  the  chattering  and  steady  state  problems  discussed  earlier. 

The  OLDTIM  AV  breakout  is  shown  in  Table  25.  The  changes  in  the  COEs  of  the 
follower  satellite  upon  which  the  OLDTIM  maneuver  is  based  are  listed  in  Table  26.  The 
OLDTIM  shows  that  as  in  the  first  simulation,  the  majority  of  the  discrete  burn  AV  is 
going  towards  the  plane  change  with  the  other  two  maneuvers  several  orders  of  magnitude 
smaller. 
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AV  (m/s) 

OLDTIM 

SMC 

Pole=-0.0006 

2.351894 

40.4172759 

Pole=-0.0012 

2.351894 

24.1352993 

Pole=-0.0018 

2.351894 

20.6293084 

Pole=-0.0024 

2.351894 

21.2748060 

Pole=-0.003 

2.351894 

23.4297615 

Table  24  Complex  Simulation  AV  Results  (Pole  Method) 


Hohmann  AV 
Inc/LAN  AV 
AP  AV 

0.00055705  (m/s) 
2.26933  (m/s) 
0.0820091  (m/s) 

Total 

2.351894  (m/s) 

Table  25  Complex  Simulation  OLDTIM  AV 

The  settling  time  results  are  listed  in  Tables  27  and  28.  The  criteria  for  settling  time 
is  that  the  position  error  for  each  direction  in  the  RCO  frame  falls  and  remains  within 
+/-  10  meters  after  which  a  stationkeeping  controller  will  take  over  and  reduce  the  error 
to  required  levels.  As  with  the  AV  costs,  the  linear  and  nonlinear  gain  controller  settling 
times  are  virtually  identically.  To  see  the  trends,  the  data  is  graphed  in  Figures  34  and 
35  for  AV  and  Figures  36  and  37  for  settling  time.  Similar  to  the  simple  p  parameter 
change,  the  gain  techniques  appear  to  asymptotically  approach  the  OLDTIM  maneuver 
as  control  usage  is  tightened  and  settling  time  increases.  The  sliding  mode  method  has 
a  local  minimum  for  both  AV*  and  settling  time  at  about  P=-0.0018  on  the  real  axis. 

All  methods  are  compared  against  each  other  and  the  OLDTIM  maneuver  in  Figure  38. 
Note  that  the  closer  the  curve  comes  to  the  origin,  the  better  the  overall  performance  of 
the  method.  Once  a  specific  method  is  determined,  a  designer  can  find  the  point  at  which 
the  tradeoff  between  AV  and  settling  time  is  acceptable  based  on  mission  parameters 
(on  board  fuel,  acceptable  off  line  times,  frequency  of  reconfigurations,  etc.).  For  the 
complex  reconfiguration,  the  gain  methods  worked  significantly  better  than  the  sliding 
mode  method.  Between  the  three  gain  methods,  there  was  no  clear  advantage  to  using  the 
nonlinear  techniques. 
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Initial 

Final 

A 

Semi-Major  Axis  (km) 

10000 

10000 

0.0017649 

Eccentricity 

0.0000499951 

0.000149848 

0.0000998524 

Inclination  (deg) 

10.0023 

9.99514 

0.00721041 

Argument  of  Perigee  (deg) 

304.987 

319.916 

14.9285 

Longitude  of  the  Ascending  Node  (deg) 

-0.00945911 

-0.120565 

0.111106 

Table  26  Complex  Simulation  Follower  COEs  Changes 


Settling  Time  (min) 

OLDTIM 

LQR 

LQR  With  LF 

SDRE 

R  =  10y 

82.933415 

21.25 

21.25 

21.25 

R  =  1010 

82.933415 

37.50 

37.50 

R  =  1011 

82.933415 

86.50 

86.50 

R  =  1012 

82.933415 

152.50 

152.50 

R  =  1013 

82.933415 

407.75 

407.25 

Table  27  Complex  Simulation  Settling  Time  Results  (Control  Weight  Methods) 


Settling  Time  (min) 

OLDTIM 

SMC 

Pole=-0.0006 

82.933415 

166.00 

Pole=-0.0012 

82.933415 

107.75 

Pole=-0.0018 

82.933415 

99.00 

Pole=-0.0024 

82.933415 

107.00 

Pole=-0.003 

82.933415 

123.00 

Table  28  Complex  Simulation  Settling  Time  Results  (Pole  Method) 


Gain  Methods 


Figure  34  Complex  Simulation  AV  (Control  Weight  Methods) 
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Pole  Methods 


Figure  35  Complex  Simulation  AV  (Pole  Method) 


Gain  Methods 


Figure  36  Complex  Simulation  Settling  Time  (Control  Weight  Methods) 
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Pole  Methods 


Figure  37  Complex  Simulation  Settling  Time  (Pole  Method) 


Method  Comparison 


Figure  38  Complex  Simulation  Method  Comparison 
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4- 5  Simple  9  Reconfiguration 

The  final  relative  orbit  reconfiguration  changed  only  the  parameter  9.  The  9  param¬ 
eter  defines  the  angle  that  the  follower  satellite’s  position  vector  makes  with  C  at  t  =  0. 
Changing  9  does  not  change  the  size,  shape,  or  orientation  of  the  relative  orbit  only  the 
follower  satellites’  starting  position  within  the  relative  orbit.  If  the  formation  included 
multiple  follower  satellites,  changing  9  would  change  that  follower’s  position  relative  to  the 
other  satellites  in  the  formation.  The  choice  of  the  initial  and  final  p  are  small  enough  that 
the  linear  Hill’s  equations  are  a  very  close  approximation  of  the  nonlinear  CW  equations. 
From  the  error  analysis  done  in  Section  2.2,  we  expect  no  more  than  a  1  meter  error  in 
the  R  and  O  directions  and  10  meters  in  the  C  direction  for  a  10,000  km  leader  orbit. 
The  initial  and  final  relative  parameters  are  given  in  Table  29.  The  simulation  uses  a  15 
second  time  step  with  the  leader  in  a  circular  orbit  with  classical  orbital  elements  (COEs) 
listed  in  Table  30.  The  AV  results  for  controllers  in  which  control  weight  was  varied  are 


Initial 

Final 

P  (km) 

0.5 

0.5 

a  (km) 

0 

0 

0  (deg) 

30 

60 

b  (km) 

0 

0 

m  (slope) 

1 

1 

n  (slope) 

0 

0 

Table  29  9  Simulation  Initial  and  Final  Parameters 


Semi-Major  Axis 

10000  km 

Eccentricity 

0 

Inclination 

10  deg 

Argument  of  Perigee 

0  deg 

Longitude  of  the  Ascending  Node 

0  deg 

Initial  True  Anomaly 

10  deg 

Table  30  9  Simulation  Leader  COEs 


listed  in  Table  31.  The  results  once  again  show  only  the  slightest  of  performance  increases 
by  using  the  nonlinear  gain  techniques  over  LQR  (on  the  order  of  10-6).  The  AV  results 
for  sliding  mode  control  (which  varied  pole  placement)  are  listed  in  Table  32.  The  slid¬ 
ing  mode  results  are  an  order  of  magnitude  higher  than  the  gain  method  controllers  with 
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AV  (m/s) 

OLDTIM 

LQR 

LQR  With  LF 

R  =  109 

R  =  1010 

R  =  1011 

R  =  1012 

R  =  1013 

0.326750 

0.326750 

0.326750 

0.326750 

0.326750 

3.6062904 

2.1439774 

1.3546992 

0.8592822 

0.5046080 

3.6062857 

2.1439539 

1.3546386 

0.8591845 

0.5044617 

3.6062818 

2.1439696 

1.3547005 

0.8592881 

0.5045870 

Table  31  6  Simulation  AV  Results  (Control  Weight  Methods) 


AV  (m/s) 

OLDTIM 

SMC 

Pole=-0.0006 

0.326750 

22.9683408 

Pole=-0.0012 

0.326750 

11.6012523 

Pole=-0.0018 

0.326750 

7.8489131 

Pole=-0.0024 

0.326750 

6.3261655 

Pole=-0.003 

0.326750 

5.2428135 

Table  32  6  Simulation  AV  Results  (Pole  Method) 

no  decrease  in  settling  time.  The  chatter  and  steady  state  problems  associated  with  the 
sliding  mode  controller  drive  up  total  AV  costs  once  again. 

Each  of  the  OLDTIM  maneuvers  is  shown  in  Table  33.  The  changes  in  the  COEs  of 
the  follower  satellite  upon  which  the  OLDTIM  maneuver  is  based  are  listed  in  Table  34. 

Hohmann  AV'  0.00000576885  (m/s) 

Inc/LAN  0.163404  (m/s) 

AP  Ay _ 0.16334  (m/s) 

Total  0.326750  (m/s) 

Table  33  9  Simulation  OLDTIM  Ay 

In  this  case  the  OLDTIM  break  out  shows  a  nearly  even  split  in  control  energy  between 
the  plane  change  (inclination  and  LAN)  and  orientation  change  (AP). 

The  settling  time  results  are  listed  in  Tables  35  and  36.  The  criteria  for  settling 
time  is  that  the  position  error  for  each  direction  in  the  RCO  frame  falls  and  remains 
within  +/-  10  meters.  It  is  assumed  that  at  this  point  a  stationkeeping  controller  will 
take  over  and  reduce  the  error  to  required  levels.  As  with  previous  examples,  the  linear 
and  nonlinear  gain  controller  settling  times  are  virtually  identically.  Finally,  the  data  is 
presented  graphically  in  Figures  39  and  40  for  Ay  and  Figures  41  and  42  for  settling 
time.  Similar  to  the  p  parameter  change,  the  gain  techniques  appear  to  asymptotically 
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Initial 

Final 

A 

Semi-Major  Axis  (km) 

10000 

10000 

0.0000182746 

Eccentricity 

0.0000499945 

0.0000499943 

2.56005  xlO"10 

Inclination  (deg) 

10.0027 

10.0018 

0.000850248 

Argument  of  Perigee  (deg) 

289.986 

319.977 

29.9911 

Longitude  of  the  Ascending  Node  (deg) 

-0.00563948 

-0.0126345 

0.00699506 

Table  34  9  Simulation  Follower  COEs 


Settling  Time  (min) 

OLDTIM 

LQR 

LQR  With  LF 

SDRE 

82.933394 

15.75 

15.75 

15.75 

82.933394 

28.75 

28.75 

28.75 

1 

82.933394 

48.75 

48.75 

48.75 

82.933394 

83.25 

83.25 

83.25 

82.933394 

180.25 

180.25 

180.25 

Table  35  9  Simulation  Settling  Time  Results  (Control  Weight  Methods) 


Settling  Time  (min) 

OLDTIM 

SMC 

Pole=-0.0006 

82.933394 

91.50 

Pole=-0.0012 

82.933394 

46.00 

Pole=-0.0018 

82.933394 

32.00 

Pole=-0.0024 

82.933394 

26.00 

Pole=-0.003 

82.933394 

22.50 

Table  36  9  Simulation  Settling  Time  Results  (Pole  Method) 
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approach  the  OLDTIM  maneuver  as  control  usage  is  tightened  and  settling  time  increases. 
The  sliding  mode  method  appears  to  approach  asymptotically  as  the  poles  on  the  real  axis 
move  farther  from  the  origin  but  a  much  slower  rate  than  the  gain  methods.  Figure 


Gain  Methods 


Figure  39  6  Simulation  AV  (Control  Weight  Methods) 

Pole  Methods 


Figure  40  6  Simulation  AV  (Pole  Method) 


43  shows  each  method  graphed  with  AV  versus  settling  time.  Note  that  the  closer  the 
curve  comes  to  the  origin,  the  better  the  overall  performance  of  the  method.  Using  this 
criteria  to  compare  methods,  it  is  easy  to  see  that  for  the  9  change  simulation  LQR,  LQR 
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Gain  Methods 


R  {Control  Weight) 

Figure  41  6  Simulation  Settling  Time  (Control  Weight  Methods) 


Pole  Methods 


Figure  42  6  Simulation  Settling  Time  (Pole  Method) 
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Method  Comparison 


Figure  43  6  Simulation  Method  Comparison 


with  linearizing  feedback,  and  SDRE  work  significantly  better  than  sliding  mode  control 
(which  once  again  must  deal  with  chattering  and  steady  state  problems).  Much  like  the 
first  simulation,  each  of  these  three  methods  produces  virtually  identical  results,  thus  there 
is  no  advantage  to  using  nonlinear  techniques  over  linear  ones. 

4.6  The  Effect  of  J2  on  Reconfiguration 

The  results  presented  assume  the  ideal  orbit  model;  the  Earth  is  a  point  source  and 
no  other  forms  of  perturbations  exist  (drag,  third  body  effects,  solar  wind,  etc.).  There  are 
two  ways  to  approach  the  perturbations  problem  in  the  relative  dynamics  frame.  The  first 
is  to  have  both  the  leader  and  follower  correct  for  perturbations  while  still  maintaining 
the  proper  relative  alignment.  The  second  is  to  accept  drift  due  to  perturbations  and  have 
the  follower  satellite  only  correct  for  the  difference  in  the  perturbations  between  the  two 
satellites.  With  either,  the  equations  of  motion  used  to  derive  the  controller  must  include 
the  perturbations  model.  Effects  caused  by  the  Earth  being  an  oblate  sphere  as  opposed  to 
a  point  mass  are  the  strongest  of  the  perturbations  encountered  by  operational  satellites. 
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Its  derivation  can  be  found  in  Appendix  B  and  in  the  IJK  frame  is 


15z2x  3x 


15  z3  _  9 z 
~dT  d? 


(119) 


where  Re  is  the  radius  of  the  Earth,  J2  is  a  constant  coefficient  found  through  observation 
of  satellite  orbits  and  is  equal  to  0.0010826  (unitless),  and  the  coordinates  of  the  satellite 
axe  x,  y,  and  z  measured  in  the  inertial  IJK  frame  ^|dj  =  d  =  \/x'2  +  y2  +  z2 j .  Thus  the 
equation  of  motion  for  a  satellite  in  inertial  space  with  only  J2  is 


+  vc  (120) 

These  modified  equations  of  motion  are  not  incorporated  into  the  development  of  the 
controllers  used  in  this  thesis  (Chapter  III),  however  simulations  were  run  with  J2  included 
into  the  inertial  propagation  of  both  leader  and  follower  to  show  its  effect.  Two  cases  are 
presented  here;  both  are  LQR  controllers  at  two  different  control  weights  ( R  =  109  and 
R  =  1013)  for  a  simple  p  change  reconfiguration  (p  =  0.5  — ►  1.5).  Both  controllers  were 
able  to  deal  with  J2  to  varying  degrees  of  success.  The  R  =  109  brought  the  error  within 
10  meters  and  the  R  =  1013  to  within  15  meters.  The  J2  perturbation  creates  a  drift 
in  the  argument  of  perigee  (AP)  and  longitude  of  the  ascending  node  (LAN)  (5).  The 
controllers  were  able  to  transform  the  AP  drift  to  a  periodic  oscillation  but  did  little  to 
correct  the  LAN  drift.  As  with  making  linear  assumptions  for  the  CW  equations,  accepting 
this  drift  term  is  an  engineering  tradeoff  for  reducing  the  complexity  of  the  system.  Given 
the  relatively  short  period  of  time  that  reconfigurations  occur  over  and  the  fact  that  only 
the  difference  in  the  J2  effect  between  the  leader  and  follower  has  to  be  incorporated  into 
the  controller,  there  is  not  an  enormous  need  to  include  J2  into  the  relative  model.  Not 
surprisingly,  the  J2  effect  prevents  the  controller  from  dampening  the  error  to  zero  and 
creates  a  periodic  error  with  an  amplitude  of  approximately  1.2  meters  for  R  =  109  and 
25  meters  for  R  =  1013.  The  small  abnormalities  can  best  be  seen  by  looking  at  the  COE 


d  = 


—pd  (  pJ2-Rg 

~ id3  ^  2 


15z2x 

d7 

15z2; 


3a; 
d * 


Lozpj  _  AQ 


15z3 

IT 


9  z 
d * 
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time  history  in  Figures  44  and  45  Note  the  scale  for  each  graph  and  that  the  oscillations 


LQR  R=1e+009 


Time  (days) 


0.05  0.1  0.15 

Time  (days) 


Figure  44  COE  Time  History  with  J2  (LQR  R  =  109) 


while  present  are  not  very  large.  Both  controllers  still  do  the  job  of  controlling  to  the 
commanded  relative  parameters  as  seen  in  Figures  46  and  47 
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V.  Conclusions  and  Recommendations 


Several  conclusions  can  be  made  from  the  data  produced  by  this  thesis.  The  first  is 
that  continuous  feedback  control  will  burn  considerably  more  fuel  than  the  discrete  burn 
method  described  in  Section  III.  Although  it  must  be  stated  again  that  the  discrete 
method  used  in  this  thesis  is  not  an  optimal  method,  the  comparison  of  it  to  the  continuous 
feedback  A  Vs  seems  to  indicate  that  there  are  more  optimal  ways  to  reconfigure  relative 
orbits.  The  feedback  controllers  do  offer  some  advantages;  namely  the  ability  to  control 
relative  position  error  to  a  leader  or  another  follower  within  the  formation.  The  biggest 
problem  with  continuous  feedback  controllers  is  that  the  controller  output  is  proportional 
to  the  error  without  regard  for  the  best  place  in  the  orbit  to  do  a  burn.  For  example,  to 
change  the  size  of  an  orbit  via  the  Hohmann  method  (an  optimal  method  for  same  plane 
circular  to  circular  orbits),  thrusting  occurs  at  apogee  and  perigee.  A  continuous  feedback 
controller  would  continue  to  thrust  throughout  the  transfer  orbit  until  the  desired  size 
change  was  attained,  burning  much  more  fuel  in  the  process.  Truly  optimal  controllers 
will  need  to  take  into  account  not  only  the  magnitude  of  the  error  to  be  corrected  but 
also  the  best  time  within  an  orbit  to  make  that  correction.  Further,  if  continuous  methods 
are  used  for  small  relative  orbits,  nonlinear  control  techniques  currently  offer  no  clear 
advantage  to  standard  linear  controllers,  and  may  in  fact  produce  inferior  results.  This 
comes  as  no  surprise  given  that  with  smaller  relative  orbits,  the  linear  Hill’s  equations 
are  a  very  good  approximation  of  the  nonlinear  CW  equations.  The  only  time  nonlinear 
controllers  prove  superior  is  when  the  relative  orbit  is  large  enough  to  create  large  errors 
between  the  linear  and  nonlinear  EOMs.  As  p  becomes  larger  (creating  a  larger  relative 
orbit),  LQR  with  LF  becomes  increasingly  more  efficient  as  compared  to  the  linear  LQR 
technique.  Since  the  command  signal  is  based  on  the  linear  EOMs,  this  simulation,  as 
implemented,  is  ideally  suited  for  small  relative  orbit  reconfigurations.  For  those  small 
orbit  reconfigurations,  LQR  (a  linear  method)  is  as  good  as  or  better  than  the  nonlinear 
methods  tested  here  as  well  as  significantly  simpler  to  implement.  Of  course,  this  thesis  did 
not  exhaust  the  numerous  nonlinear  control  techniques  applicable  to  this  problem.  These 
controllers  should  be  explored  before  nonlinear  techniques  are  completely  discarded.  Also, 
in  order  to  model  more  realistic  conditions,  measurement  noise  should  be  added  and/or 


70 


an  estimator  used  instead  of  assuming  full  state  feedback.  Finally,  there  is  great  research 
potential  in  exploring  the  nonlinear  EOMs  in  order  to  develop  command  signals  that  are 
valid  for  larger  relative  orbits.  Perhaps  a  Keplerian  type  formulation  of  the  nonlinear 
EOMs  based  on  an  elliptical  orbit.  If  nothing  else,  techniques  need  to  be  developed  to 
choose  appropriate  initial  conditions  such  that  the  secular  drift  term  in  the  C  direction  is 
nulled  out. 
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Appendix  A.  The  Two  Body  Problem 

To  find  the  acceleration  of  a  satellite  relative  to  a  much  larger  primary  body  due  to  the 
force  of  gravity,  consider  Figure  48  where  d\  and  d<i  are  position  vectors  in  inertial  space. 


/V 


Figure  48  Two  Body  Problem 


By  Newton’s  second  law  (19) 


y;  F  =  mass  *  acc 

—4  -4  -4 

Fg  +  fc  +  fp  =  mass  *  acc 


(121) 

(122) 


where  Fg  is  the  force  due  to  gravity,  fc  is  a  vector  of  control  forces,  and  fp  is  a  vector  of 
perturbative  forces.  Gravitational  forces  follow  an  inverse  square  law  (19).  Using  Figure 
48 

Mass  i  *  Mass2 


\Fg\  °c 


I  d< 


(123) 


To  make  this  an  equality  the  universal  constant  of  gravitation  G  is  used  (G  =  6.672  x  10  11 

G  *  Mass\  *  Mass2 


\Fg 


I  d? 


(124) 
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Subbing  into  Equation  122  for  both  masses 


G  *  Massi  *  Mass2 


L  l<*l2 

G  *  Mass\  *  Mass2 

.  w 


U  +  /ci  +  fpi  —  Massi  *  d\ 


(“ U )  +  fc2  +  fp2  —  MaSS2  *  c?2 


where  U  is  a  unit  vector  parallel  to  d  and  gives  the  scaler  force  of  gravity  a  direction. 
Making  the  following  substitution 

U  =  i  (126) 

\d\ 


yields 


G  *  Massi  *  Mass2  d  ,  y  ,  ?  ..  y 

- HE -  *  7=f  +  /ci  +  fpi  =  Massi  *  di 

\d\ 2  J  |d| 

G  *  Massi  *  Mass2  d  -  -  y 

- -=5T -  *~y+fc2  +  Jp2  =  Mass2  *  d2 

I «  J  l«l 


Combining  terms  and  dividing  through  by  the  masses  to  isolate  the  acceleration 


G  *  Mass2  y  /ci  ,  fpi 

- — -  %  Ct  ~ r -  ~j~ - 

|rf|3  J  Mass  i  Massi 

-G*  Massi]  y  fc2  ,  fp 2 

- =5 -  *  d  +  — - f-  — - 

|d|3  Mass2  Massi 


Since  d  =  c?2  —  di,  taking  the  derivative  twice  yields 


d  =  c?2  —  di 


Subbing  Equations  128  into  the  above  equation 


y  — G  *  Massi  y  fc2  fp 2 

d  =  - =7- -  *«+— - b  -77 - 

|d  3  J  Mass2  Mass2 


G^Massi  (  (130) 

|  d  3  Mass  1  Massi 


Combining  terms 


j*=  —G{Massi  +  Mass2)d  _  /ci _ fpi  fc 2  /P2  /131n 

|Jj3  Massi  Massi  Mass2  Mass2 
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If  first  body  is  the  Earth  and  the  second  body  a  satellite  orbiting  the  Earth  then  Mass\  = 
Mass  Earth)  MaSS2  —  Mass  Sat 5  /cl  ”  f EarthCcmtroU  fc2  =  f SatControl)  fpl  =  fEarthPert)  <Uld 
fp2  =  fSatPert 


G ( Af as S Earth  "b  Masssat )d 

w 


+ 


f SatControl 

Masssat 


fSatPert 

Masssat 


f EarthControl 

Mass  Earth 


fEarthPert 

Mass  Earth 


(132) 


Since  MassEarth  >  Masssat  we  can  say  that  Mass^^  +  Mass^at  «  MassEarth ,  further, 
if  the  Earth  is  being  used  as  the  inertial  reference  frame  for  the  satellite  then  f  EarthControl  = 
0  and  fEarthPert  =  0* 


"j*  — G  *  Mass  Earth  *  d  .  f SatControl  .  fSatPert 

d  =  - — -  +  ~1T7 - 1- 


\d\ 3 


Masssat  Masssat 


(133) 


If  we  define  the  constant  fi 


fjt  =  G*  MassEarth  =  6.672  x  KT11^  x  5.974236  x  10 24kg 

H  =  3.98601  x  1014^  =  398601^  (134) 


and  define  the  specific  forces 


►  f SatControl 

Masssat 

V  —  jjatPert 

p  Masssat 


(135) 

(136) 


Then  the  equation  of  motion  for  a  satellite  about  the  Earth  with  control  input  of  vc  and 
a  perturbative  acceleration  of  vp  is 

d  =  H—  +  vc  +  Up  (137) 

I  d  3 
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Appendix  B.  The  Ji  Perturbation 

Because  gravity  is  a  conservative  force,  it  can  be  derived  from  the  gradient  of  a  scalar 
potential  function  (20).  Looking  at  Equation  137  without  any  control  inputs 

{m) 

Expanding  in  the  inertial  reference  frame  shown  in  Figure  48  and  noting  that  |d|  =  d  — 
y/x 2  +  y2  +  «‘2  where  x,  y,  and  z  are  the  position  of  the  satellite  in  the  inertial  frame  IJK 


d  =  - <==■ - s+vpi  I  +  - — - T  +  VpJ  J 

\x2  +  y2  +  z2]  2  J  [[x2+y2  +  z2]  2 


Noting  that 

dip 


dx  \MI  J  L 

—  fix 

[x2  +  y2  +  z2]  § 

where  B  is  a  potential  function  such  that 


+  - 5-  +  VpK  K 

[[X2 +  y2  +  Z2]?  J 


z2  +  rf  +  z2*  (0)  -  /i[l/2(\/a;2  +  y2  +  ^)-1/2]2x 
x2  +  y2  +  z2 

—  fix 


+  VpJ 


+  Vpi 


Similarly 


9_  (ft 
dy  \d 


VB  =  vp 

(141) 

=  21,  +v 

(142) 

\x2  +  y2  +  z2]  2 

—uz 

=  .  „  „  „.3  +VpK 

(143) 

[x2  +  y2  +  z2]  2 

£(n +b)  =  — 

dz\d  )  [x2  +  92  +  22]1  PK 

Therefore  the  inertial  acceleration  is  the  gradient  of  the  potential  function  (j  +  B 


d  =  V  q  +  B 
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The  B  term  comes  about  by  modeling  the  Earth  not  as  a  point  mass  but  as  an  oblate  body 
with  nonhomogeneous  mass  distribution  (Figure  49)  and  is  the  sum  on  the  infinite  series 


Figure  49  Oblate  Earth  Coordinate  Frame 


{OO  /  \n  n  ^  /  Jl  \n  ) 

Y  [-J-J  JnPn(sin<f>)  +  Y  (C'nm  cos  (p  +  Snm sin  cp)Pnm(s'm <j>)  | 


<p  =  m\  +  U)ete 


(145) 

(146) 


where  Re  is  the  mean  equatorial  radius  of  the  Earth,  0  is  the  geocentric  latitude  of  the 
satellite  (measured  from  the  equator),  A  is  the  geographical  longitude  (measured  from  the 
prime  meridian),  d  is  the  satellite  position  vector  magnitude,  u>e  is  the  rotation  rate  of  the 
Earth,  te  is  the  time  since  the  I  direction  lined  up  with  the  Greenwich  meridian,  Jn  is  the 
zonal  harmonic  coefficients  of  order  0,  Pn  is  a  Legendre  polynomial  of  degree  n  and  order 
0,  Pnm  is  a  Legendre  polynomial  of  degree  n  and  order  m ,  Cnm  is  the  tesseral  harmonic 
coefficient  for  and  Snm  is  the  sectorical  harmonic  coefficents  for  n  =  m. 

Measurements  of  the  zonal,  tesseral  and  sectorial  coefficients  (Jn,  Cnrn.  and  Snm) 
reveal  that  the  J2  term  is  at  least  400  times  larger  than  the  next  most  significant  term. 
Thus  for  applications  such  as  satellite  reconfigurations  which  occur  over  relatively  short 
time  periods  all  higher  terms  can  be  ignored.  Using  this  assumption,  Equation  145  reduces 
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B= -it [(f)  J”p”<sinH  =  -i (f )  w™*) 

where  J2  =  0.0010826  and  the  2nd  Legendre  polynomial  has  the  form  (12) 


P2pO  =  -(3X2-l) 


Using  this  Legendre  polynomial  in  Equation  147 


b=zt{t)2 


rearranging  terms  in  B 


By  geometry  in  Figure  49 


B=^(9i[3<8in«2-10 


sin  (j)  —  - 
a 


which  means 


_  J2R2  (  fj,  f  z2 

B~^r [w 


1  = 


/. iJ2R2  f  —Sz2  +  d2 


Now  that  B  is  in  cartesian  coordinates,  vp  is  found  by  taking  the  gradient  of  B 


*P  =  VB= 


d5(2df)-(-3z2+<i2)5(i4§ 
2Rl  d5(2rf2i)-(-3z2+d2)5d4* 

2  d m 
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Simplifying 
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Appendix  C.  Clohessy-  Wiltshire  Derivation 

Assume  the  relationship  between  the  leader  and  follower  satellite  shown  in  Figure  50  and 
that  the  leader  satellite  is  in  a  circular  orbit  about  the  Earth  (Section  1.1).  In  the  RCO 


Figure  50  Clohessy-Wiltshire  Setup 
frame  (Figure  1)  the  position  vectors  are: 

(155) 

(156) 

(157) 


L  —  k0R  +  0C*  +  0  O 
P  =  rR  +  cC  +  oO 
M  =  L  +  P  =  {k0  +  r)R  +  cC  +  06 


where  k0  is  the  orbit  radius  of  the  leader  satellite.  To  take  the  inertial  derivative  of  a 
vector: 


_RCO 


d(X) 

_|_  jyRCO—wrt— Inertial  x 


(158) 


where  N  is  the  angular  velocity  of  the  RCO  frame  with  respect  to  the  inertial  frame,  in 
this  case  NRCO-WTt-Inertml  =  QR+OC +uiO  where  u  is  the  angular  frequency  of  the  leader 
satellite’s  orbit  and  is  found  via  Equation  3.  Taking  the  inertial  derivative  of  Equation  157 
yields  the  velocity: 

M  =RCO  M  +  (wO)  x  M  (159) 


where  RCOM  is  the  derivative  of  M  in  the  RCO  frame. 


rcoM  =  rR  +  cC  +  60 


(160) 
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R  C 

6 

— ucR 

(ud)  x  M  = 

0  0 

UJ 

= 

uj{k0  +  t)C 

kQ  +  r  c 

o 

o  6 

Adding  these  together  yields: 

M  =  (r  —  uc)R  +  [c  +  ui(k0  +  r)]C  +  oO 


(161) 


(162) 


Taking  the  inertial  derivative  of  Equation  162  yields  the  acceleration: 

k  =RCO  k  +  (u6)  x  k 

where  RCOM  is  the  derivative  of  M  in  the  RCO  frame. 

rcoM  =  (r  —  ujc)R  +  (c  +  ujf)C  +  '60 


R  C 

6 

[— uc  —  u2(k0  +  r)]R 

(ud)  xM  = 

0  0 

UJ 

— 

(uf  -  <jJ2c)C 

r  —  u>c  c  +  u(k0  +  r ) 

o 

o  6 

Adding  these  together  yields: 

M  =  [f  -  2 uc  -  u>2(k0  +  r)]E  +  (c  +  2ur  -  u2c)C  +  60 


(163) 


(164) 

(165) 


(166) 


Prom  Appendix  A  we  know  that  acceleration  due  to  gravity  and  control  (no  perturbations) 
is: 

..TUT 

(167) 


A  -uM  _ 
M  =  — +  vc 


\M\ 

where  /i  is  the  gravitational  constant  and  v*c  is  the  control  acceleration.  By  Kepler  s  third 
law  (29) 

uj  =  4  /■/?  M  =  u2kl  (168) 
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Expanding  the  magnitude  of  the  position  vector 


|M|3  = 


\j{k0  +  r) 


2  +  c2  +  o2 


(K  +  r)2  +  c2  +  o2 


(169) 


and  using  these  equalities  in  Equation  167,  we  find  that 

~  —u>2k%[(k0  +  r)R  +  cC  +  oO]  _ 

M  = - — - - - 3 - -  +  vc 

[(k0  +  r)2  +  c2  +  o2]a 


(170) 


Equating  166  and  170 

[r  —  2uc  —  ui2{k0  +  r)].R 
(c  +  2  ujr  —  uj2c)C 
60 


[(k0  +  r )2  +  c2  +  o2]  2 


{k0  +  r)R 
cC 
oO 


+ 


VcrR 

VccC 

vc00 


(171) 


Equating  scalar  terms  and  simplifying  yields  the  CW  equations: 


f  —  2o;c  —  u2{k0  +  r ) 


1  - 


[(fco+r)2+c2+02]l 


-vCr=0 


c  +  2  uir  —  uj2c 
d  +  u>2o 


1- 


[(fco+r)2+c2+02]2  J 

UCo  =  0 


“  VCc  =  0 


L[(fco+02+c2+°2]5  J 


(172) 
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Appendix  D.  Linear  Clohessy- Wiltshire  Derivation  (Hill’s  Equations) 
Starting  with  the  CW  Equations  from  Appendix  C 


r  —  2  uic  —  u2(k0  +  r) 


1  - 


£ 


[(fco+r)2+e2+o2]2  J 


c  +  2uir  —  uj2c 
d  +  ui2o 


1  - 


kl 


[(feo+r)2+c2+o2]5  J 


-  VCr  =  0 

vCc  =  0 


jg_ 


L[(feo+r)2+c2+o2]^  J 


-  Vco  =  0 


(173) 


The  nonlinear  part  of  these  equations  is 


NL  =  [(&o  +  r)2  +  c2  +  o2]  2  =  [k2  +  2  kQr  +  r2  +  c2  +  o2]  2 


(174) 


Factoring  out  an  k2 


NL  = 


l2/  2r  r2  c2  o2N 

M1  +  F  +  F  +  fc2  +  A-2 
\  ft'0  '*'0  1^0  , 


" 2  3  /  2r  r2  c2  o2  \  2 

=  ko[1  +  T  +  T2  +  u2+T2) 

V  k°  ko  kb  kB  J 


(175) 


Using  binomial  expansion  (13),  this  becomes 


NL  =  kz0 


1  +  l(£  +  FI  +  *?)+H°'T 


(176) 


Neglecting  higher  order  terms  (H.O.T)  and  assuming  that  r,  c,  and  o  are  appropriately 
small  compared  to  the  radius  of  the  leader  orbit  (Chapter  1.1)  such  that 


r2  c 2  o 2 

fc2  ~fc2  ~fc2  ~° 
Ko 


(177) 


With  this  assumption,  the  nonlinear  term  becomes 


(178) 


Subbing  into  Equations  173 


r  —  2ujc  — u>2(k0  +  r)  1  —  ^ 


+3  r 


-Vcr  =0 
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(179) 


C+2ujf-  U)2C  [l  -  5^]  -  VCc  =  0 

°  +  u2°  [fc^3f]  ~vco=0 

With  some  algebra 

r-2uc-  u2(k0  +  r )  ~vCr=0 

•c  +  2ur-u2c\j^\  -vCc  =  0  (180) 

d  +  uj2°  [*&f]  ~vc0  =  0 

Assuming  that  r  is  appropriately  small  such  that 


k0  +  r  s=s  k0 

(181) 

k0  +  3r  «  kQ 

(182) 

yeilds 

r  -  2 uc  -  u2k0  -  Vcr  =  0 

c  +  2 uir  -  u2  -  vCc  =  0 

(183) 

Assuming  that  r  and 

o  +  u2o  —  vCo  =  0 

c  are  appropriately  small  such  that 

cr 

Y0~ 

(184) 

Yeilds  the  linearized  1 

version  of  the  CW  equations 

r  —  2  uc  —  3  u2r  —  vCr  =  0 

c  +  2  wr  —  vCc  =  0 

(185) 

d  +  u2o  -vCo=0 
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Appendix  E.  Parametrization  of  Hill’s  Equations 

Starting  with  the  homogeneous  Hill’s  equations  from  Appendix  D 


r  —  2  ojc  —  3a )2r  =  0 

c  +  2a of  =  0  (186) 

o  +  a;2o  =  0 


Taking  the  Laplace  transform 

[s2il(s)  -  sr0  -  r0]  -  2u[sC{s)  -  c0]  —  3w2J ?(s)  =  0 

[s2C(s)  -  sc0  —  c0 1  +  2o;[si?(s)  —  r0\  =  0  (187) 

[s20(s)  —  so0  —  d0]  +  u>20(s )  =  0 


where  r0  and  r0  are  the  position  and  velocity  initial  conditions  (ICs)  in  the  radial  direction, 
c0  and  c0  are  the  position  and  velocity  ICs  in  the  cross  track  direction,  and  o0  and  o0  are 
the  position  and  velocity  ICs  in  the  out  of  plane  direction.  Collecting  terms  and  placing 
the  equations  in  matrix  form 


s2  —  3cj2  —2ujs 

0 

'  R(s) 

sr0  +  r0  -  2 uc0 

2  us  s 2 

0 

C(s) 

— 

sc0  +  60  +  2u >r0 

i - 

o 

o 

S2  +  U)2 

.  0 (*)  . 

so0  +  o0 

Finding  the  inverse  and  solving  for  R(s),  C(s),  0(s ) 


'  m ' 

i 

2uj 

0 

sr0  +  r'0  -  2o;c0 

S2+CJ2 

s(s2+cj2) 

C(s) 

_ 

—2u> 

s2— 3a;2 

0 

sc0  +  c0  +  2a >r0 

s(sU+UJ2) 

s2(s:l+ai2) 

.  °(s) . 

0 

0 

1 

so0  +  o0 

S2-\-(jJ2 

(188) 


(189) 
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Multiplying  this  out 


'  R(s) 

srn+r„—2<jjc„  i  2u;(sco+Co+2u/r0) 

%4  ZP  + - s^s^+o^) - 

C(s) 

— 

— 2u>(sr0+f0— 2ojc0)  .  (s2-3cj2)(sc0+c'o+2ajro) 

s(s2-fa>2)  S2(s24-CJ2) 

.  0 . 

L  S2+u4  J 

(190) 


Performing  partial  fraction  expansion 


’  R(s)  ' 

srn+fo~2(jJC0  i  2u2 Co -2scn—4sujrn  |  2co+4u;r,0  1 
s2+a;2  *  u(s2+d2)  1  U)  s 

C(s) 

= 

4sc0+4c0+8u;ro  (  -~2u2rn+2srn— 4scrta;  ,  -3c0-6cjrn  •  a;c0-2ro  1 
s2+a;2  a;(s2*fa;2)  s2  a ;  s 

.  °(s) . 

J 

(191) 


Collecting  terms 


'  i?(s)  ‘ 

^J.o_2ca^wra)s  (£)w  2cVi+4u>rv>  1  1 

s2+a;2  s2+u;2  a;  s 

C(s) 

= 

2("^')s  |  2(  °  cj  a  r°)a;  3a;  2c0+4a;r0  1  .  ucn-2r0  1 

s2+a;2  s2+a;2  2  a;  s2  a ;  s 

son  i  (l? 

s2-fa>2  '  s2+u;2 

(192) 


Taking  the  inverse  Laplace  transform 


r(t) 

j r0  —  -Cg~^-— aJ  cos(a;i  +  0)  +  *£  sin (ut  +  6)  +  2cp~^a;rp 

c(i) 

= 

2  (£)  cos  (ut  +  9)  +  2  ro]  sin(a;t  +  0)  3f  (2Cp+^rp)  + 

.  . 

Oo  cos(cj£  +  #)  +  sin(a;t  +  0) 

(193) 


where  0  is  an  initial  angle  at  t  =  0.  Since  the  ICs  and  u  do  not  vary  with  time,  two 
constants  can  be  defined 


2  c0  +  4 <jr0 
a  = - 

UJ 


u 


(194) 

(195) 


which  simplifies  Equation  193  to 


r(t )  =  (r0  —  a)  cos(ujt  +  6)  +  —  sin (u)t  +  9)  +  a 

id 


(196) 
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(197) 


c(t)  =  cos  (u)t  +  9)  +  2  (a  —  rQ )  sin(u;£  +  9)  —  +  f> 

u;  2 

o(t)  =  o0  cos  (ut  +  9)  +  —  sin(a>£  +  9)  (198) 

UJ 

For  simplicity,  r  =  r(f),  c  =  c(t),  o  —  o{t).  Isolating  the  sinusodal  terms  in  the  first  two 


equations. 


r  —  a  =  (ra  —  a)  cos(u it  +  9)  +  —  sin(wf  +  9) 

UJ 

,  3 mat  _  u  • 

- 2 - =  —  cos  (ut  +  9)  +  (a  —  r0 )  sin(u>f  +  9) 

2  UJ 


Squaring  both  sides 


[r  —  a]2  =  ( r0  —  a)2  cos 2(ut  +  9)  +  (  —  J  sm2(ut  +  9) 


[c+^fl-6]2 


+  2  (r0  —  a)—  cos(a>£  +  9)  sin(u;£  +  9) 

UJ 

=  cos2(wt  +  6)  +  (— l)2(r0  —  a)2  sin2(wf  +  9) 

+  2(— l)(r0  —  a)—  cos(wt  +  9)  sin(a>£  +  9) 

UJ 


Adding  these  equations  together 


[r  -  a]2  + 


[c  +  ^fl  -bf 


=  (r0  —  a)2  +  |  cos 2(uit  +  9)  + 


-  )  +  (r0  —  a)2  sin2(w<  +  9) 


Pulling  out  the  common  coefficient  and  using  the  identity  cos2(r)  +  sin2(r)  =  1 


r  ,  Zojat  _  k]2 

[r-a]2+[C+  24  &J 


(  \2  ,  (^°\ 
=  (r0- a)  +  — 

V  u>  J 


The  right  side  of  this  equation  has  no  time  varying  terms,  thus  define  the  constant  p 


(r„  -  af  +  (^ 
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Subsisting  and  dividing  by  p2  yields 


[r  —  a]2 
P 2 


+ 


[c+  ZfZ-b]2 

4  p2 


=  1 


(206) 


This  is  the  standard  equation  for  an  ellipse,  thus  the  follower  satellite  will  follow  an  elliptical 
path  about  the  leader  satellite  in  the  R/C  plane  (Figure  51).  Note  that  for  a  stable  orbit 


Figure  51  General  Follower  Orbit  in  the  R/C  Plane 

(i.e.  one  that  does  not  drift  over  time)  about  the  leader  satellite,  a  must  equal  zero. 
Ellipses  can  be  parameterized  in  the  general  form 

x  —  h  =  >lsin(a) 
y  —  k  =  B  cos(a) 

where  A  is  the  semi  major  or  minor  axes  in  the  x  direction,  B  is  the  semi  major  or  minor 
axes  in  the  y  direction,  and  the  center  of  the  ellipse  is  located  at  (h,k).  Placing  Equation 
206  in  this  form 

r(t )  —  a  =  p  sin(a;t  +  6)  (207) 

c(t)  +  —  b  =  2p  cos(u>t  +  0)  (208) 

z 
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and  solving  for  r(t)  and  c(t) 


r(t )  =  psin(c<;t  +  6)  +  a 


(209) 


c(t )  =  2pcos(u>t  +  9) - - — I-  b  (210) 

A 

If  we  visualize  the  elliptical  orbit  from  Figure  51  rotating  about  the  major  or  minor  axes, 
it  is  easy  to  see  it  will  trace  a  line  in  R/O  plane  (if  rotated  about  the  major  axis)  and  the 
R/C  plane  (if  rotated  about  the  minor  axis).  Using  the  equation  of  a  line  and  adding  the 
effects  of  a  rotation  about  both  major  and  minor  axes 

o(t )  =  m[r(t)  —  a]  +  n  c(t)  +  —  b  (211) 

z 

where  m  and  n  are  the  slopes  of  the  lines  formed  by  the  rotation  about  the  minor  and 
major  axes  respectively.  Subbing  Equations  207  and  208  into  the  above  equation  yields: 

o(t )  =  mpsm(u;t  +  6)  +  2npcos(uit  +  0)  (212) 

This  is  the  correct  form  expected  from  the  inverse  Laplace  transform  found  in  Equation 
198  ,  an  oscillatory  function  with  a  constant  amplitude  and  a  period  of  ut. 

Of  the  six  constants  required  to  define  a  relative  orbit  about  a  leader  satellite,  only 
three  have  been  defined  so  far  (a,  b,  and  p).  The  other  three  need  to  be  developed  in  terms 
of  initial  conditions.  Using  Equation  207  and  the  derivative  of  Equation  209  (at  t  =  0) 

rQ  —  a  —  psin(0)  (213) 

r'0  =  pucos(d)  <$■?£  =  pcos(0)  (214) 

and  inserting  them  into  Equation  212  and  it’s  derivative  (again  at  t  =  0) 

o0  =  mpsin(0)  +  2npcos(6)  (215) 

o0  =  mpucos(9 )  —  2npu>sin(0)  (216) 


87 


yields 


Oo  =  m(r0  —  a)  +  n  (2^) 
d0  =  mr'o  -  n[2u(r0  -  a)] 


In  matrix  form 


Op 

rp~  a  ^ 

m 

Op 

1 

O 

1 

to 

o 

l 

1 _ 

n 

Solving  for  m  and  n  by  multiplying  both  sides  by  the  matrix  inverse 


o 

1 

cs 

3 

1 

To 

- 1 

o° 

_ 1 

m 

u2r2-2r0w2a+ijJ2a2+ro2 

u2r2—2ruj2a+uJ2a2+fo2 

n 

i(a-r0)w 

- 1 

*o 

u;2r2 -2r0u2a+w2a2+fo2 

u)2r2—  2r0u;2a+u>2a2+fo2 

Multiplying  through  and  simplifying 


_  dpTo  -  OpUJ^a  -  r0) 
m  r'02 +u2(a-r0)2 
_  OoTpUJ  +  o0uj(a  -  rp) 
2 [r02  +  v2(a  -  r0)2] 

Finally,  to  find  9  divide  Equation  213  by  Equation  214 


rp  -  a 

ISL 
u ; 

and  solving  for  9 

Tp-a)' 

Tp 

All  parameters  can  be  solved  given  the  relative  position  and  velocity  except  for  9 , 
find  9  at  any  given  time,  divide  Equation  207  by  its  derivative  r  =  pu  cos(u )t  +  6) 

j~a  =  tan  (ut  +  9) 

LJ 


6  =  arctan 


uj( 


(217) 

(218) 

(219) 

(220) 

(221) 

(222) 

(223) 

(224) 
so  to 

(225) 
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using  the  tan  addition  identity  to  expand  the  right  side  of  the  equation  (26) 


r  —  a  _  tan(wt)  +  tan(0) 
£  1  —  tan(wt)  tan(0) 


(226) 


With  some  algebra 


u(r  —  a)  —  r  tan(wt)  =  tan(0)  [u;(r  —  a)  tan(wt)  +  r]  (227) 


and  solving  for  9 


6  =  arctan 


u>(r  —  a)  —  rtan(wt) 
w(r  —  a)  tan(wt)  +  r 


which  at  t  =  0  reduces  to  Equation  224. 


(228) 


E.l  Unit  Analysis 

Since  the  initial  and  final  relative  orbits  are  defined  in  terms  of  the  relative  parameters 
derived  above,  it  is  important  to  understand  the  proper  units  for  each  one.  The  variables 
r,  c,  and  o  are  positions  and  are  defined  in  terms  of  length.  Likewise,  r,  c,  and  6  are 
velocities  and  are  expressed  as  Finally,  to  is  the  angular  frequency  of  the  orbit  and 

is  expressed  in 


a  = 


^  +  mslen9th 
1 

time 


length 

-lift-  =  length 
time 


(229) 


b  = 


^^9th  - 

1 

time 


length 

time 


length 

=  length 

time 


(230) 


P  = 


\ 


(i length  —  length )2  + 


'  length  ' 
time 

1 

.  time  , 


=  ijlength?  +  length?  =  length 


(231) 


length  length  lrnaf  h  (  1  ) 

time  time  {time ) 

2 

|  ( length  —  length )  | 

( length  \ 
^  time  ) 

, 2  length 2 

1  time 2 

l 

(®)  +{tiL)  (len9th  ^ngth)2  | 

f  length  \ 
\  time  ) 

1 2  length 2 

'  time 2 

Oi length?  (!g^)2  -  !|Sgi 


=  unitless 

(232) 

=  unitless 

(233) 
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6  =  arctan 


[length  -  length) 

length 

time 


— — time  ~  arctan 
time 


"  length ' 
time 
length 
.  time  - 


=  degrees 


(234) 


In  the  simulations  run  for  this  thesis,  a,  6,  and  p  are  defined  in  kilometers,  9  in  degrees, 


and  m  and  n  in  terms  of  a  unitless  slope. 


90 


Appendix  F.  Inertial  to  Relative  Position  and  Velocity  Transformation 

The  transformation  matrix  from  the  inertial  frame  to  the  relative  frame  is  simply  the 
transpose  of  Equation  14.  Let  IJKL,  IJKL  and  IJKM>  IJKM  be  the  position  and  velocity 
of  the  leader  and  follower  satellites  respectively  expressed  in  the  UK  frame. 


qRC02IJK 


=  [£|C0] 


where 


ft  ijkl 
\IJKL\ 

c  =  IJKl 
0  =  RxC 


taking  the  transpose  to  get  the  transformation  from  UK  to  RCO 


qIJK2RCO 


[C 


<RC02IJK'\T 


\R\C\0]1 


(235) 


(236) 

(237) 

(238) 


(239) 


Note  that  for  orthonormal  frames  CIJK2RCO*CRCO‘2IJK  equals  the  identity  matrix.  Equa¬ 
tion  18  can  then  be  solved  for  the  relative  position 


_  qIJK2RCO^IJK m  _IJK 


(240) 


Solving  Equation  25  for  relative  velocity  {u>  is  calculated  with  Equation  3) 


r 

—  OUJ 

c 

o 

_  qIJK2RCO\IJK _IJK  _ 

ru 

0 

(241) 
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Appendix  G.  LQR  Controller  Simulation  Plots 

This  appendix  presents  more  detailed  data  from  the  simulation  discussed  in  Section  4.4 
for  an  LQR  controller.  For  convenience,  Table  37  lists  the  initial  and  final  relative  param¬ 
eters.  The  reconfiguration  is  commanded  at  0.05  days  (72  minutes)  from  the  start  of  the 
simulation.  The  control  weight  was  varied  in  order  to  produce  different  responses  for  the 
LQR  controller.  High  control  weight  translates  to  less  control  usage  and  lower  AV  with  a 
higher  settling  time  and  the  reverse  for  smaller  control  weights.  To  show  contrast,  the  low 
(R  =  109)  and  the  high  (R  =  1013)  runs  are  presented. 


Initial 

Final 

P  (km) 

0.5 

1.5 

a  (km) 

0 

0 

0  (deg) 

45 

60 

b  (km) 

0 

0 

m  (slope) 

1 

1.5 

n  (slope) 

0 

1 

Table  37  Complex  Simulation  Initial  and  Final  Parameters 

G.l  LQRR=  109 

For  the  complex  simulation,  the  LQR  controller  with  a  gain  of  R  =  109  required 
31.4307824  m/s  of  Ay  and  settled  to  within  +/-  10  meters  in  each  axis  after  21.25  minutes. 
The  motion  of  the  follower  satellite  in  the  leader  centered  RCO  frame  is  shown  in  Figure 
52  where  the  follower  starts  on  the  inner  ellipse  and  is  commanded  onto  the  outer  ellipse. 
This  change  in  relative  parameters  translates  to  a  change  in  the  inertial  COEs  (Figure  53) 
Changes  in  each  of  the  COEs  is  listed  in  Table  26.  The  control  usage  profile  necessary  to 
guide  the  follower  onto  the  desired  orbit  is  displayed  in  Figure  54.  Note  that  the  control 
inputs  are  shown  in  the  inertial  UK  frame.  The  absolute  value  of  the  control  inputs  are 
integrated  in  order  to  get  the  total  Ay  expended;  their  time  history  is  shown  in  Figure  55 
The  position  error  which  represents  half  of  the  error  signal  (the  other  half  is  the  velocity 
error)  and  upon  which  the  settling  time  is  based  is  shown  in  Figure  56.  Finally,  the  relative 
parameters  (upon  which  the  command  signal  is  based)  are  shown  as  a  function  of  time  in 
Figure  57. 
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Z  Axis  Control  Input  (m/s2)  X  Axis  Control  Input  (m/s2) 


LQR  R=1  e+009 


Time  (days) 

Figure  56  Position  Error  Time  History  (LQR  R  =  109) 
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Figure  57  Relative  Parameter  Time  History  (LQR  R  =  109) 
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G.2  LQRR=lOn 


For  the  complex  simulation,  the  LQR  controller  with  a  gain  of  R  =  1013  required 
4.2184381  m/s  of  AV  and  settled  to  within  +/- 10  meters  in  each  axis  after  407.25  minutes. 
The  motion  of  the  follower  satellite  in  the  leader  centered  RCO  frame  is  shown  in  Figure 
58  where  the  follower  starts  on  the  inner  ellipse  and  is  commanded  onto  the  outer  ellipse. 


Motion  of  Follower  Satellite  wrt  Leader  Sat  (LQR  R=1e+013) 


R  Hat  (km)  C  Hat  (km) 

Figure  58  Path  of  Follower  Satellite  in  RCO  Frame  (LQR  R  =  1013) 

This  change  in  relative  parameters  translates  to  a  change  in  the  inertial  COEs  (Figure  59) 
Changes  in  each  of  the  COEs  is  listed  in  Table  26.  The  control  usage  profile  necessary  to 
guide  the  follower  onto  the  desired  orbit  is  displayed  in  Figure  60.  Note  that  the  control 
inputs  are  shown  in  the  inertial  UK  frame.  The  absolute  value  of  the  control  inputs  are 
integrated  in  order  to  get  the  total  AV  expended;  their  time  history  is  shown  in  Figure  61 
The  position  error  which  represents  half  of  the  error  signal  (the  other  half  is  the  velocity 
error)  and  upon  which  the  settling  time  is  based  is  shown  in  Figure  62.  Finally,  the  relative 
parameters  (upon  which  the  command  signal  is  based)  are  shown  as  a  function  of  time  in 
Figure  63. 
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Figure  61  AV  Time  History  (LQR  R  =  1013) 


Figure  62  Position  Error  Time  History  (LQR  R  =  1013) 
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Figure  63  Relative  Parameter  Time  History  (LQR  R  =  1013) 
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Appendix  H.  LQR  with  LF  Controller  Simulation  Plots 

This  appendix  presents  more  detailed  data  from  the  simulation  discussed  in  Section  4.4 
for  an  LQR  with  linearizing  feedback  (LF)  controller.  For  convenience,  Table  38  lists 
the  initial  and  final  relative  parameters.  The  reconfiguration  is  commanded  at  0.05  days 
(72  minutes)  from  the  start  of  the  simulation.  The  control  weight  was  varied  in  order  to 
produce  different  responses  for  the  LQR  with  LF  controller.  High  control  weight  translates 
to  less  control  usage  and  lower  AF  with  a  higher  settling  time  and  the  reverse  for  smaller 
control  weights.  To  show  contrast,  the  low  ( R  =  109)  and  the  high  (R  =  1013)  runs  are 
presented. 
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Table  38  Complex  Simulation  Initial  and  Final  Parameters 


H.l  LQR  with  LF  R  =  109 

For  the  complex  simulation,  the  LQR  with  LF  controller  with  a  gain  of  R  =  109 
required  31.4307616  m/s  of  AV  and  settled  to  within  +/-  10  meters  in  each  axis  after 
21.25  minutes.  The  motion  of  the  follower  satellite  in  the  leader  centered  RCO  frame  is 
shown  in  Figure  64  where  the  follower  starts  on  the  inner  ellipse  and  is  commanded  onto 
the  outer  ellipse.  This  change  in  relative  parameters  translates  to  a  change  in  the  inertial 
COEs  (Figure  65)  Changes  in  each  of  the  COEs  is  listed  in  Table  26.  The  control  usage 
profile  necessary  to  guide  the  follower  onto  the  desired  orbit  is  displayed  in  Figure  66.  Note 
that  the  control  inputs  are  shown  in  the  inertial  UK  frame.  The  absolute  value  of  the 
control  inputs  are  integrated  in  order  to  get  the  total  AV  expended;  their  time  history  is 
shown  in  Figure  67  The  position  error  which  represents  half  of  the  error  signal  (the  other 
half  is  the  velocity  error)  and  upon  which  the  settling  time  is  based  is  shown  in  Figure  68. 
Finally,  the  relative  parameters  (upon  which  the  command  signal  is  based)  are  shown  as  a 
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Motion  of  Follower  Satellite  wrt  Leader  Sat  (LQR  WLF  R=1e+009) 
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Figure  64  Path  of  Follower  Satellite  in  RCO  Frame  (LQRWLF  R  =  109) 
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Figure  65  Time  History  of  the  Follower  Satellite  COEs  (LQRWLF  R  =  109) 
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LQR  WLF  R=1e+009 
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Figure  68  Position  Error  Time  History  (LQRWLF  R  =  109) 
function  of  time  in  Figure  69. 
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LQR  WLF  R=1  e+009 
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Figure  69  Relative  Parameter  Time  History  (LQRWLF  R  =  109) 

H.2  LQR  with  LFR=  1013 

For  the  complex  simulation,  the  LQR  with  LF  controller  with  a  gain  of  R  =  1013 
required  4.2181246  m/s  of  AV  and  settled  to  within  +/-  10  meters  in  each  axis  after 
407.75  minutes.  The  motion  of  the  follower  satellite  in  the  leader  centered  RCO  frame  is 
shown  in  Figure  70  where  the  follower  starts  on  the  inner  ellipse  and  is  commanded  onto 
the  outer  ellipse.  This  change  in  relative  parameters  translates  to  a  change  in  the  inertial 
COEs  (Figure  71)  Changes  in  each  of  the  COEs  is  listed  in  Table  26.  The  control  usage 
profile  necessary  to  guide  the  follower  onto  the  desired  orbit  is  displayed  in  Figure  72.  Note 
that  the  control  inputs  axe  shown  in  the  inertial  UK  frame.  The  absolute  value  of  the 
control  inputs  are  integrated  in  order  to  get  the  total  AV  expended;  their  time  history  is 
shown  in  Figure  73  The  position  error  which  represents  half  of  the  error  signal  (the  other 
half  is  the  velocity  error)  and  upon  which  the  settling  time  is  based  is  shown  in  Figure  74. 
Finally,  the  relative  parameters  (upon  which  the  command  signal  is  based)  are  shown  as  a 
function  of  time  in  Figure  75. 
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Motion  of  Follower  Satellite  wrt  Leader  Sat  (LQR  WLF  R=1e+013) 


R  Hat  (km)  C  Hat  (km) 

Figure  70  Path  of  Follower  Satellite  in  RCO  Frame  (LQRWLF  R  =  1013) 
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Figure  71  Time  History  of  the  Follower  Satellite  COEs  (LQRWLF  R  =  1013) 
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Figure  72  Control  Usage  (LQRWLF  R  =  1013) 
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Figure  73  AU  Time  History  (LQRWLF  R  =  1013) 
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Figure  74  Position  Error  Time  History  (LQRWLF  R  =  1013) 
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Figure  75  Relative  Parameter  Time  History  (LQRWLF  R  =  1013) 
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Appendix  I.  SDRE  Controller  Simulation  Plots 

This  appendix  presents  more  detailed  data  from  the  simulation  discussed  in  Section  4.4  for 
the  State  Dependent  Riccati  Equations  (SDRE)  controller.  For  convenience,  Table  39  lists 
the  initial  and  final  relative  parameters.  The  reconfiguration  is  commanded  at  0.05  days 
(72  minutes)  from  the  start  of  the  simulation.  The  control  weight  was  varied  in  order  to 
produce  different  responses  for  the  SDRE  controller.  High  control  weight  translates  to  less 
control  usage  and  lower  AV  with  a  higher  settling  time  and  the  reverse  for  smaller  control 
weights.  To  show  contrast,  the  low  (R  =  109)  and  the  high  (R  =  1013)  runs  are  presented. 
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Table  39  Complex  Simulation  Initial  and  Final  Parameters 

1.1  SDRER  =  109 

For  the  complex  simulation,  the  SDRE  controller  with  a  gain  of  R  =  109  required 
31.4307437  m/s  of  AV  and  settled  to  within  +/-  10  meters  in  each  axis  after  21.25  minutes. 
The  motion  of  the  follower  satellite  in  the  leader  centered  RCO  frame  is  shown  in  Figure 
76  where  the  follower  starts  on  the  inner  ellipse  and  is  commanded  onto  the  outer  ellipse. 
This  change  in  relative  parameters  translates  to  a  change  in  the  inertial  COEs  (Figure  77) 
Changes  in  each  of  the  COEs  is  listed  in  Table  26.  The  control  usage  profile  necessary  to 
guide  the  follower  onto  the  desired  orbit  is  displayed  in  Figure  78.  Note  that  the  control 
inputs  are  shown  in  the  inertial  UK  frame.  The  absolute  value  of  the  control  inputs  are 
integrated  in  order  to  get  the  total  AV  expended;  their  time  history  is  shown  in  Figure  79 
The  position  error  which  represents  half  of  the  error  signal  (the  other  half  is  the  velocity 
error)  and  upon  which  the  settling  time  is  based  is  shown  in  Figure  80.  Finally,  the  relative 
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Motion  of  Follower  Satellite  wrt  Leader  Sat  (SDRE  R=1e+009) 
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Figure  76  Path  of  Follower  Satellite  in  RCO  Frame  (SDRE  R  =  109) 
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Figure  77  Time  History  of  the  Follower  Satellite  COEs  (SDRE  R  =  109) 
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Z  Axis  Control  Input  (m/s2) 
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Figure  80  Position  Error  Time  History  (SDRE  R  =  109) 

parameters  (upon  which  the  command  signal  is  based)  are  shown  as  a  function  of  time  in 
Figure  81. 
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Figure  81  Relative  Parameter  Time  History  (SDRE  R  =  109) 


1.2  SDRER  =  1013 

For  the  complex  simulation,  the  SDRE  controller  with  a  gain  of  R  =  1013  required 
4.2185121  m/s  of  AV  and  settled  to  within  +/- 10  meters  in  each  axis  after  407.25  minutes. 
The  motion  of  the  follower  satellite  in  the  leader  centered  RCO  frame  is  shown  in  Figure 
82  where  the  follower  starts  on  the  inner  ellipse  and  is  commanded  onto  the  outer  ellipse. 
This  change  in  relative  parameters  translates  to  a  change  in  the  inertial  COEs  (Figure  83) 
Changes  in  each  of  the  COEs  is  listed  in  Table  26.  The  control  usage  profile  necessary  to 
guide  the  follower  onto  the  desired  orbit  is  displayed  in  Figure  84.  Note  that  the  control 
inputs  are  shown  in  the  inertial  UK  frame.  The  absolute  value  of  the  control  inputs  are 
integrated  in  order  to  get  the  total  AV  expended;  their  time  history  is  shown  in  Figure  85 
The  position  error  which  represents  half  of  the  error  signal  (the  other  half  is  the  velocity 
error)  and  upon  which  the  settling  time  is  based  is  shown  in  Figure  86.  Finally,  the  relative 
parameters  (upon  which  the  command  signal  is  based)  are  shown  as  a  function  of  time  in 
Figure  87. 
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Figure  82  Path  of  Follower  Satellite  in  RCO  Frame  (SDRE  R  =  1013) 
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Figure  83  Time  History  of  the  Follower  Satellite  COEs  (SDRE  R  =  1013) 
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Figure  84  Control  Usage  (SDRE  R  =  1013) 
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Figure  85  AV  Time  History  (SDRE  R  =  1013) 
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Figure  86  Position  Error 
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Appendix  J.  SMC  Controller  Simulation  Plots 

This  appendix  presents  more  detailed  data  from  the  simulation  discussed  in  Section  4.4 
for  the  sliding  mode  controller.  For  convenience,  Table  40  lists  the  initial  and  final  relative 
parameters.  The  reconfiguration  is  commanded  at  0.05  days  (72  minutes)  from  the  start 
of  the  simulation.  The  pole  position  was  varied  in  order  to  produce  different  responses  for 
the  sliding  mode  controller.  Poles  farther  from  the  ju>  axis  translate  to  more  control  usage 
and  higher  AV  with  a  lower  settling  time  and  the  reverse  for  smaller  control  weights.  To 
show  contrast,  the  P  =  —0.0006  and  P  =  —0.0018  runs  are  presented. 
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Table  40  Complex  Simulation  Initial  and  Final  Parameters 

J.l  SMC  P  =  -0.0006 

For  the  complex  simulation,  the  SMC  controller  with  a  gain  of  P  =  —0.0006  required 
40.4172759  m/s  of  AV  and  settled  to  within  +/-  10  meters  in  each  axis  after  166  minutes. 
The  motion  of  the  follower  satellite  in  the  leader  centered  RCO  frame  is  shown  in  Figure 
88  where  the  follower  starts  on  the  inner  ellipse  and  is  commanded  onto  the  outer  ellipse. 
This  change  in  relative  parameters  translates  to  a  change  in  the  inertial  COEs  (Figure  89) 
Changes  in  each  of  the  COEs  is  listed  in  Table  26.  The  control  usage  profile  necessary  to 
guide  the  follower  onto  the  desired  orbit  is  displayed  in  Figure  90.  Note  that  the  control 
inputs  are  shown  in  the  inertial  UK  frame.  The  absolute  value  of  the  control  inputs  are 
integrated  in  order  to  get  the  total  AV  expended;  their  time  history  is  shown  in  Figure  91 
The  position  error  which  represents  half  of  the  error  signal  (the  other  half  is  the  velocity 
error)  and  upon  which  the  settling  time  is  based  is  shown  in  Figure  92.  Finally,  the  relative 
parameters  (upon  which  the  command  signal  is  based)  are  shown  as  a  function  of  time  in 
Figure  93. 
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Figure  89  Time  History  of  the  Follower  Satellite  COEs  (SMC  P 
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Figure  92  Position  Error  Time  History  (SMC  P  =  — OJ 
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Figure  93  Relative  Parameter  Time  History  (SMC  P  =  — 0.( 


51 


J.2  SMC  P  =  -0.0018 


For  the  complex  simulation,  the  SDRE  controller  with  a  gain  of  P  =  —0.0018  required 
20.6293084  m/s  of  AV  and  settled  to  within  +/-  10  meters  in  each  axis  after  99  minutes. 
The  motion  of  the  follower  satellite  in  the  leader  centered  RCO  frame  is  shown  in  Figure 
94  where  the  follower  starts  on  the  inner  ellipse  and  is  commanded  onto  the  outer  ellipse. 


Motion  of  Follower  Satellite  wrt  Leader  Sat  (NS  SMC  Pole=-0.0018) 
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Figure  94  Path  of  Follower  Satellite  in  RCO  Frame  (SMC  P  =  —0.0018) 

This  change  in  relative  parameters  translates  to  a  change  in  the  inertial  COEs  (Figure  95) 
Changes  in  each  of  the  COEs  is  listed  in  Table  26.  The  control  usage  profile  necessary  to 
guide  the  follower  onto  the  desired  orbit  is  displayed  in  Figure  96.  Note  that  the  control 
inputs  are  shown  in  the  inertial  UK  frame.  The  absolute  value  of  the  control  inputs  are 
integrated  in  order  to  get  the  total  AV  expended;  their  time  history  is  shown  in  Figure  97 
The  position  error  which  represents  half  of  the  error  signal  (the  other  half  is  the  velocity 
error)  and  upon  which  the  settling  time  is  based  is  shown  in  Figure  98.  Finally,  the  relative 
parameters  (upon  which  the  command  signal  is  based)  are  shown  as  a  function  of  time  in 
Figure  99. 


120 


Axis  Control  Input  (m/s2)  X  Axis  Control  Input  (m/s2)  LAN(degs)  Inc(degs) 


NS  SMC  Pole=-0.0018 


10000.5 

^  10000 

I 

<  9999.5 
W  9999 
9998.5 

0  0.1  0.2  0.3  0.4 

10.01 

10 

9.99 

9.98 

( 

0 

-0.05 
-0.1 
-0.15 
-0.2 

0  0.1  0.2  0.3  0.4 

Time  (days) 


Time  (days) 


Time  (days) 


Time  History  of  the  Follower  Satellite  COEs  (SMC  P  = 


Time  (days) 

Figure  96  Control  Usage  (SMC  P  =  —0.0018) 
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Appendix  K.  LQR  Simulink  Diagram 
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Appendix  L.  LQR  With  Linearizing  Feedback  Simulink  Diagram 


Figure  101  LQR  With  Linearizing  Feedback  Simulink  Diagram 
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Appendix  M.  State  Dependent  Riccati  Equations  Simulink  Diagram 
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Appendix  N.  Sliding  Mode  Control  Simulink  Diagram 
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